Estadística No Paramétrica

Universidade de Oriente, Venezuela

1 Pruebas que Usan Datos de una Muestra Aleatoria Simple

En este núcleo temático se describen de manera teórica algunas pruebas no paramétricas que usan una muestra aleatoria simple para hacer inferencias acerca de los parámetros de la población de la cual fue tomada. Dentro de estas pruebas se encuentran:

1. La Prueba Binomial (Sección 1.1)

2. La Prueba del Signo (Sección 1.2)

3. La Prueba de Bondad de Ajuste Ji-cuadrado (Sección 1.4)

4. La Prueba de Rangos Asignados de Wilcoxon (Sección 1.3)

5. La Prueba de Bondad de Ajuste Kolmogorov-Smirnov (Sección 1.5)

6. La Prueba de Rachas o Carreras de Walf-Wolfowitz (Sección 1.6)

Una vez descrita estas pruebas se procede, mediante ejemplos, a la implementación en R de las mismas.

El siguiente bloque de código permite instalar y cargar los paquetes usados para implementar estas pruebas en R.

packages <- c(
  "rmarkdown", "bookdown", "bookdownplus", "magrittr",
  "kableExtra", "BSDA", "DT", "devtools", "tidyverse",
  "htmltools", "htmlwidgets", "utf8", "stringr", "MASS",
  "Hmisc", "kolmim", "hrbrthemes", "plotly", "reshape2",
  "randtests", "snpar", "knitcitations", "exactRankTests",
  "stats", "exact2x2", "NSM3", "gmodels", "RVAideMemoire",
  "MultNonParam"
)
package.check <- lapply(packages, FUN = function(x) {
  if (!require(x, character.only = TRUE)) {
    install.packages(x, dependencies = TRUE)
    library(x, character.only = TRUE)
  }
})

1.1 Prueba Binomial

La prueba binomial es una prueba no paramétrica utilizada para contrastar hipótesis referente a la proporción de individuos de una poblacional que poseen una característica de interés. La aplicación de este test requiere que la variable objeto de estudio sea categórica dicotómica o se haya dicotomizado a través de alguna transformación.

1.1.1 Racionalización de la Prueba Binomial

En muchos problemas de la vida real el estudio se centra en analizar variables aleatorias de tipo categóricas dicotómicas (o que se han dicotomizado a través de alguna transformación). Suponga que una de las categorías se codifica como éxito de acuerdo si el valor observado en la variable se corresponde con una característica de interés, de lo contrario se codifica como fracaso. Si se toma una muestra aleatoria cualquiera de tamaño \(n\) de esta variable, tal que:

\[ X_{i}=\begin{cases} 1, & \textit{si el i-ésimo valor observado en la muestra es un éxito}\\ 0, & \textit{si el i-ésimo valor observado en la muestra es un fracaso} \end{cases}, \]

con \(i=1, 2, \dotsc,n\). Suponiendo que para cada valor observado en la muestra, la probabilidad de que éste sea un éxito es constante e igual a \(p\), es decir, \((P(X_i=1)=p)\), lo que implica que la probabilidad de observar un fracaso es \(q=1-p\). Visto desde otro punto de vista, \(p\) puede ser interpretado como la proporción de individuos en la población, de la cual se extrajo la muestra, que poseen la característica de interés; lo cual se justifica si la muestra se toma con reemplazo, si la población es infinita o si la población es finita tal que la aproximación de la distribución hipergeométrica a la distribución binomial sea adecuada. En base a lo anterior, se define el estadístico de prueba

\[ \begin{equation} X=\sum_{i=1}^{n}X_{i}, \tag{1.1} \end{equation} \]

el cual representa el número de éxitos o individuos que poseen la característica de interés en la muestra. Bajo los supuestos establecidos anteriormente, este estadístico tiene una distribución binomial con parámetros \(n\) y \(p\), es decir, \(X\sim B\left ( n,p \right )\). En tal sentido, la probabilidad de observar de manera exacta \(x\) éxitos en una muestra viene dada por la función de probabilidad binomial, la cual se muestra a continuación:

\[ \begin{equation} P\left ( X=x \right )=p_{B}\left ( x,\ n,\ p \right )=\binom{n}{x}p^{x}q^{n-x}, \ x= 0, 1, 2, \dotsc, n, \ n>0 \ y \ 0<p<1. \tag{1.2} \end{equation} \]

Mientras que la probabilidad de obtener a lo más \(x\) éxitos se puede calcular a través de la función de distribución acumulada de probabilidad binomial, la cual se muestra en la siguiente ecuación:

\[ \begin{equation} P\left ( X\leq x \right )=F_{B}\left ( x,\ n,\ p \right )=\sum_{i=0}^{x}\binom{n}{i}p^{i}q^{n-i}, \ x= 0, 1, 2, \dotsc, n, \ n>0 \ y \ 0<p<1. \tag{1.3} \end{equation} \]

Por otro lado, si se quiere contrastar la hipótesis nula de que la proporción de individuos en la muestra que poseen la característica de interés es \(p_{0}\) \(\left(H_0: p=p_0\right)\) contra la hipótesis alternativa unilateral superior de que la proporción de individuos en la población que poseen esta característica es mayor que la establecida en \(H_{0}\), es decir, \(H_1:p>p_{0}\). En base a la información de una muestra, una estimación de \(p\) daría luces al respecto. Se sabe que \(\widehat{P}=\frac{X}{n}\) es el estimador de máxima verosimilitud de \(p\), por lo tanto si \(x\) es el valor de \(X\) observado en la muestra, \(\widehat{p}=\frac{x}{n}\) es una estimación de \(p\). Ahora, si se supone que la hipótesis nula es cierta \(\widehat{p}\) debe ser próximo a \(p_{0}\), de lo contrario habría indicios de que esta hipótesis es falsa. Una forma de medir esta distancia, en términos relativos, entre \(\widehat{p}\) y \(p_{0}\) es a través de la probabilidad de obtener un valor \(\widehat{P}\) tan grande o más extremo que \(\widehat{p}\) (\(p_{valor}\)), es decir, \(P\left ( \widehat{P}\geq \widehat{p} \right )=p_{valor}\). Pero como la distribución de \(\widehat{P}\) es desconocida, se usa la distribución de \(X\) para calcular el \(p_{valor}\). Dado que \(X\) se distribuye binomialmente, bajo la hipótesis nula, con parámetros \(n\) y \(p=p_{0}\), el \(p_{valor}\) se calcula de la siguiente manera:

\[ p_{valor}=P\left ( \widehat{P} \geq \widehat{p}\right )=P\left ( \frac{X}{n} \geq \frac{x}{n}\right )=P\left ( X\geq x \right )=1-P\left ( X\leq x-1 \right ). \]

Dado que, bajo la hipótesis nula, \(X\sim B\left ( n,p = p_{0} \right )\), de la ecuación (1.3), el \(p_{valor}\), viene determinado por:

\[ \begin{equation} p_{valor}=1-F_{B}\left ( x-1, \ n, \ p=p_{0} \right )=1-\sum_{i=0}^{x-1}\binom{n}{i}p^{i}q^{n-i}. \tag{1.4} \end{equation} \]

En caso de que la hipótesis alternativa sea unilateral inferior, es decir, \(H_{1}:p<p_{0}\), el \(p_{valor}\) representa la probabilidad de encontrar un valor de \(\widehat{P}\) tan pequeño o más pequeño \(\widehat{p}\). O sea,

\[ p_{valor}=P\left ( \widehat{P} \leq \widehat{p}\right )=P\left ( \frac{X}{n} \leq \frac{x}{n}\right )=P\left ( X\leq x \right )=P\left ( X\leq x \right ). \]

Lo que se puede reescribir, según la ecuación (1.3), de la siguiente manera,

\[ \begin{equation} p_{valor}=F_{B}\left (x, \ n, \ p=p_{0} \right )=\sum_{i=0}^{x}\binom{n}{i}p^{i}q^{n-i}. \tag{1.5} \end{equation} \]

Si el \(p_{valor}\leq\alpha\) se acepta \(H_{1}\), ya sea que esta haya sido planteado en forma unilateral superior o unilateral inferior.

Para una hipótesis alternativa de dos colas, \(H_{1}=p \neq p_0\) , el \(p_{valor}\) viene dado por:

\[ \begin{equation} p_{valor}= \begin{cases} 2\left [ 1-F_{B}\left ( x-1,\ n,\ p_{0} \right ) \right ], \ \textit{si} \ x> np_{0} \\ 2F_{B}\left ( x,\ n, \ p_{0} \right ), \ \textit{si} \ x\leq np_{0} \end{cases}. \tag{1.6} \end{equation} \]

Note que el cálculo del \(p_{valor}\) depende de la evaluación de la función de distribución binomial \(\left ( F_{B} \left ( x,n,p_{0} \right )\right )\), la cual se encuentra ampliamente tabulada para diferentes valores de \(n\) y \(p\) (generalmente \(n\leq25\) y \(p\leq1/2\)). En la actualidad estas tablas han quedado en desuso dado que existen una gran cantidad de softwares informáticos que permiten evaluar esta función sin prácticamente restricciones sobre estos parámetros. En este sentido, la distribución base de R proporciona las funciones:

dbinom(x, size, prob, log = FALSE)
pbinom(q, size, prob, lower.tail = TRUE, log.p = FALSE)
qbinom(p, size, prob, lower.tail = TRUE, log.p = FALSE)
rbinom(n, size, prob)

La función dbinom permite evaluar la función de masa de probabilidad binomial definida por la ecuación (1.2), pbinom calcula la función de distribución binomial descrita por la ecuación (1.3), qbinom determina los cuantiles binomiales y rbinom genera números pseudoaleatorios binomiales. Los parámetros de estas funciones se describen en la tabla 1.1.

knitr::kable(
  data.frame(
    stringsAsFactors = TRUE,
    Parámetros = c(
      "x, q", "p", "n", "size", "prob", "log, log.p", "lower.tail"
    ),
    Descripción = c(
      "Vector de cuantiles.",
      "Vector de Probabilidades.",
      "Número de observaciones. Si length(n) > 1, se considera que la longitud es el número requerido.",
      "Número de ensayos (cero o más).",
      "Probabilidad de éxito en cada ensayo.",
      "Lógico; si es TRUE, las probabilidades p se dan como log(p).",
      "Lógico; si es TRUE (por defecto), las probabilidades son $P[X \\leq x]$, de lo contrario, $P[X > x]$."
    )
  ),
  booktabs = TRUE,
  caption = "Descripción de los parámetros de funciones de R relacionadas con la distribución binomial",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE,
    htmltable_class = "lightable-hover"
  ) %>%
  kable_classic_2()

Tabla 1.1: Descripción de los parámetros de funciones de R relacionadas con la distribución binomial

Parámetros	Descripción
x, q	Vector de cuantiles.
p	Vector de Probabilidades.
n	Número de observaciones. Si length(n) > 1, se considera que la longitud es el número requerido.
size	Número de ensayos (cero o más).
prob	Probabilidad de éxito en cada ensayo.
log, log.p	Lógico; si es TRUE, las probabilidades p se dan como log(p).
lower.tail	Lógico; si es TRUE (por defecto), las probabilidades son \(P[X \leq x]\), de lo contrario, \(P[X > x]\).

También, la distribución base de R proporciona la función binom.test(x, n, p = 0.5, alternative = c("two.sided", "less", "greater"), conf.level = 0.95); la cual ejecuta la prueba binomial exacta. La descripción de los parámetros de esta función se muestran en la tabla 1.2.

knitr::kable(
  data.frame(
    stringsAsFactors = FALSE,
    Parámetros = c("x", "n", "p", "alternativa", "conf.level"),
    Descripción = c(
      "Número de éxitos, o un vector de longitud 2 que indique el número de éxitos y fracasos, respectivamente.",
      "Número de ensayos; se ignora si x tiene longitud 2.",
      "La probabilidad de éxito especificada en la hipótesis nula",
      "Indica la hipótesis alternativa y debe ser: \"two.sided\",
      \"greater\" o \"less\". Puede especificar sólo la letra inicial.",
      "Nivel de confianza para el intervalo de confianza establecido."
    )
  ),
  booktabs = TRUE,
  caption = "Descripción de los parámetros de la función binom.test de R",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2()

Tabla 1.2: Descripción de los parámetros de la función binom.test de R

Parámetros	Descripción
x	Número de éxitos, o un vector de longitud 2 que indique el número de éxitos y fracasos, respectivamente.
n	Número de ensayos; se ignora si x tiene longitud 2.
p	La probabilidad de éxito especificada en la hipótesis nula
alternativa	Indica la hipótesis alternativa y debe ser: “two.sided”, “greater” o “less”. Puede especificar sólo la letra inicial.
conf.level	Nivel de confianza para el intervalo de confianza establecido.

Como se dijo anteriormente, el cálculo del \(p_{valor}\) depende de la evaluación de la función de distribución acumulada binomial dada en la ecuación (1.3), lo que implica un gran volumen de cálculo cuando el tamaño de la muestra es grande. Pero, dado que \(X\) es una suma de \(n\) variables aleatorias independientes de Bernoulli, por el teorema del límite central, la distribución de \(X\) se aproxima a la normal en la medida que \(n\) tiende a infinito. Esta convergencia se acelera en la medida que \(p\) tienda \(\frac{1}{2}\). Al respecto, Canavos (1988, pág. 142), establece como criterio de una adecuada aproximación si \(np>5\) cuando \(p≤\frac{1}{2}\) o si \(nq>5\) cuando \(p>\frac{1}{2}\). Como \(X\) se distribuye binomialmente, bajo la hipótesis nula, su media viene dada por \(\mu_{X}=np_{0}\) y su varianza por \(\sigma_{X}=np_{0}q_{0}\). En consecuencia, la variable aleatoria \(X\simeq N\left ( \mu_{X}=np_{0}, \sigma_{X}=\sqrt{np_{0}q_{0}} \right )\). Por lo que:

\[ \begin{equation} Z=\frac{X-np_{0}}{\sqrt{np_{0}q_{0}}}, \tag{1.7} \end{equation} \]

tiene distribución aproximadamente normal estándar, es decir, \(Z\simeq N\left(\mu_{Z}=0, \sigma_{Z}=1\right)\), para un \(n\) suficientemente grande. Note que lo que se intenta es aproximar probabilidades de una variable discreta a partir de una variable aleatoria continua, por lo que la aproximación puede mejorarse aplicando la corrección por continuidad de Yate. En este caso el estadístico \(Z\) se redefine de la siguiente manera,

\[ \begin{equation} Z=\frac{\left(X\pm 0,5 \right)-np_{0}}{\sqrt{np_{0}q_{0}}}, \tag{1.8} \end{equation} \]

donde \(X + 0,5\) se usa cuando \(X < np\) y \(X - 0,5\) se usa cuando\(X > np\).

En el caso de una hipótesis alternativa unilateral superior \(\left(H_1:p>p_{0}\right)\), la aproximación del \(p_{valor}\) se obtiene de la siguiente manera:

\[ p_{valor}=P\left ( X\geq x \right )= 1-P\left ( X\leq x-1 \right )\cong 1-P\left ( Z\leq \frac{x-0,5-np_{0}}{\sqrt{np_{0}q_{0}}} \right ). \]

En resumen,

\[ \begin{equation} p_{valor}\cong 1-F_{N}\left ( Z= \frac{x-0,5-np_{0}}{\sqrt{np_{0}q_{0}}}, \ \mu_{Z}=0, \ \sigma_{Z}=1 \right ). \tag{1.9} \end{equation} \]

Ahora, cuando la hipótesis alternativa es unilateral inferior \(\left(H_1: p<p_0\right)\), la aproximación del \(p_{valor}\) se obtiene, usando la corrección de Yate, como sigue:

\[ p_{valor}=P\left(X \leq x\right) \cong P\left(X≤x+0,5\right). \]

Es decir,

\[ \begin{equation} p_{valor}\cong F_{N}\left ( Z= \frac{x+0,5-np_{0}}{\sqrt{np_{0}q_{0}}}, \ \mu_{Z}=0, \ \sigma_{Z}=1 \right ). \tag{1.10} \end{equation} \]

Una aproximación del \(p_{valor}\), cuando la prueba de hipótesis es bilateral \(\left(H_1: p\neq p_0\right)\), se obtiene con la siguiente ecuación

\[ \begin{equation} p_{valor}\cong \begin{cases} 2\left [ 1-F_{N}\left ( Z=\frac{x-0,5-np_{0}}{\sqrt{np_{0}q_{0}}}, \ \mu_{Z}=0, \ \sigma_{Z}=1\right ) \right ] , \ \textit{si} \ x>np_{0} \\ 2 F_{N}\left ( Z= \frac{x+0,5-np_{0}}{\sqrt{np_{0}q_{0}}} ,\ \mu_{Z}=0, \ \sigma_{Z}=1\right ), \ \textit{si} \ x \leq np_{0} \end{cases}. \tag{1.11} \end{equation} \]

Como se puede observar, el calculo del \(p_{valor}\) aproximado depende de la función de distribución normal estándar, como se observa en las ecuaciones (1.9), (1.10) y (1.11). La distribución base de R proporciona esta y otras funciones relacionadas con la distribución normal, tales como:

dnorm(x, mean = 0, sd = 1, log = FALSE)
pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)
rnorm(n, mean = 0, sd = 1)

La tabla 1.3 describe los parámetros de estas funciones

knitr::kable(
  data.frame(
    stringsAsFactors = FALSE,
    Parámetros = c(
      "x, q", "p", "n", "mean",
      "sd", "log, log.p", "lower.tail"
    ),
    Descripción = c(
      "Vector de cuantiles.",
      "Vector de Probabilidades.",
      "Número de observaciones. Si length(n) > 1, se considera que la longitud es el número requerido.", "Vector de medias.",
      "Vector de desviaciones estándar.",
      "Lógico; si es TRUE, las probabilidades p se dan como log(p).",
      "Lógico; si es TRUE (por defecto), las probabilidades son $P[X \\leq x]$, de lo contrario, $P[X > x]$."
    )
  ),
  booktabs = TRUE,
  caption = "Descripción de los parámetros de funciones de R relacionadas con la distribución normal",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2()

Tabla 1.3: Descripción de los parámetros de funciones de R relacionadas con la distribución normal

Parámetros	Descripción
x, q	Vector de cuantiles.
p	Vector de Probabilidades.
n	Número de observaciones. Si length(n) > 1, se considera que la longitud es el número requerido.
mean	Vector de medias.
sd	Vector de desviaciones estándar.
log, log.p	Lógico; si es TRUE, las probabilidades p se dan como log(p).
lower.tail	Lógico; si es TRUE (por defecto), las probabilidades son \(P[X \leq x]\), de lo contrario, \(P[X > x]\).

La distribución base de R incluye la función prop.test(x, n, p = NULL, alternative = c("two.sided", "less", "greater"), conf.level = 0.95, correct = TRUE), la cual permite calcular el \(p_{valor}\) aproximado de la prueba binomial para muestra grandes. Los parámetros de esta función se especifican en la tabla 1.4.

knitr::kable(
  data.frame(
    stringsAsFactors = FALSE,
    Parámetros = c("x", "n", "p", "alternative", "conf.level", "correct"),
    Descripción = c(
      "Un vector de conteos de éxitos, una tabla unidimensional con dos entradas, o una tabla bidimensional (o matriz) con dos columnas, dando los conteos de éxitos y fracasos, respectivamente.",
      "Un vector de conteos de ensayos; ignorado si x es una matriz o una tabla.",
      "Un vector de probabilidades de éxito. La longitud de p debe ser igual al número de grupos especificados por x, y sus elementos deben ser mayores que 0 y menores que 1.",
      "Una cadena de caracteres que especifique la hipótesis alternativa, que puede ser: \"two.sided\" (por defecto), \"greater\" o \"less\". Puede especificar sólo la letra inicial. Sólo se utiliza para comprobar que una sola proporción es igual a un valor dado, o que dos proporciones son iguales; en caso contrario, se ignora.",
      "Nivel de confianza del intervalo de confianza establecido. Debe ser un solo número entre 0 y 1. Sólo se usa cuando se comprueba que una sola proporción es igual a un valor dado, o que dos proporciones son iguales; en caso contrario se ignora.",
      "Una indicación lógica de si la corrección de continuidad de Yates debe aplicarse siempre que sea posible."
    )
  ),
  booktabs = TRUE,
  caption = "Descripción de los parámetros de la función prop.test de R",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 1.4: Descripción de los parámetros de la función prop.test de R

Parámetros	Descripción
x	Un vector de conteos de éxitos, una tabla unidimensional con dos entradas, o una tabla bidimensional (o matriz) con dos columnas, dando los conteos de éxitos y fracasos, respectivamente.
n	Un vector de conteos de ensayos; ignorado si x es una matriz o una tabla.
p	Un vector de probabilidades de éxito. La longitud de p debe ser igual al número de grupos especificados por x, y sus elementos deben ser mayores que 0 y menores que 1.
alternative	Una cadena de caracteres que especifique la hipótesis alternativa, que puede ser: “two.sided” (por defecto), “greater” o “less”. Puede especificar sólo la letra inicial. Sólo se utiliza para comprobar que una sola proporción es igual a un valor dado, o que dos proporciones son iguales; en caso contrario, se ignora.
conf.level	Nivel de confianza del intervalo de confianza establecido. Debe ser un solo número entre 0 y 1. Sólo se usa cuando se comprueba que una sola proporción es igual a un valor dado, o que dos proporciones son iguales; en caso contrario se ignora.
correct	Una indicación lógica de si la corrección de continuidad de Yates debe aplicarse siempre que sea posible.

1.1.2 Implementación en R de la Prueba Binomial

Ejemplo 1.1 (Prueba Binomial) En un estudio de los efectos del estrés, un investigador enseño a 18 estudiantes dos métodos diferentes para hacer el mismo nudo. La mitad de los sujetos (seleccionados aleatoriamente) aprendieron primero el método A y la otra mitad aprendió primero el método B. Posteriormente –a medianoche y después de un examen final de cuatro horas de duración–, a cada sujeto se le pidió que hiciera el nudo. La predicción fue que el estrés induciría regresión, esto es, que los sujetos regresarían al primer método aprendido para hacer el nudo. Cada sujeto fue categorizado conforme a si usó el primero o el segundo método aprendido de hacer nudos cuando se le pedía que hiciera el nudo bajo estrés. Resultando que 16 de los 18 utilizó el primer método aprendido para hacer el nudo. Este ejemplo es tomado de Siegel (1995, pág. 64).

En este ejemplo, la hipótesis alternativa que se quiere contrastar es que el estrés induce regresión en los estudiantes, es decir, que bajo estrés, la proporción de estudiantes que realiza el primer método aprendido para hacer el nudo es mayor que los que usan el segundo método aprendido \(\left(H_1:p>\frac{1}{2}\right)\). Dicho de otra manera, bajo estrés, la proporción de estudiantes que deciden realizar el primer nudo aprendido es mayor que 0,5.

El siguiente código muestra la evaluación de la ecuación (1.4) usando la función pbinom, la cual evalúa la función de distribución acumulada binomial con parámetros \(n=18\) y \(p=0,5\), en \(X=16-1=15\); es decir, \(F\left(X=15, n=18,p=0,5\right)\).

x <- 16
n <- 18
p0 <- 0.5
p.valor <- 1 - pbinom(q = x - 1, size = n, prob = p0)
paste("El p.valor = ", p.valor)

#> [1] "El p.valor =  0.0006561279296875"

Como se puede ver en la salida generada por el código anterior, el \(p_{valor}= 0,0007\) es menor que el nivel de significancia establecido (\(\alpha = 0,05\)), por lo que se concluye que el estrés produce regresión en los estudiantes, es decir, bajo estrés, la mayoría de los estudiantes deciden realizar el primer nudo aprendido.

El resultado anterior, también se puede obtener a través de la función binom.test, cuyos parámetros se describen en la tabla 1.2, a través del siguiente código:

x <- 16
n <- 18
p0 <- 0.5
(ejemplo1 <- binom.test(
  x = x, n = n, p = p0, alternative = "greater",
  conf.level = 0.95
))

#> 
#>  Exact binomial test
#> 
#> data:  x and n
#> number of successes = 16, number of trials = 18, p-value = 0.0006561
#> alternative hypothesis: true probability of success is greater than 0.5
#> 95 percent confidence interval:
#>  0.6897373 1.0000000
#> sample estimates:
#> probability of success 
#>              0.8888889

Como se nota en la salida, el \(p_{valor}\) obtenido por este método es igual al obtenido por el procedimiento anterior.

Por otro lado, la aproximación asintótica del \(p_{valor}\) indicado por la ecuación (1.9), se puede obtener con la ejecución del siguiente código.

x <- 16
n <- 18
p0 <- 0.5
z <- (x - 0.5 - n * p0) / sqrt(n * p0 * (1 - p0))
(p.valor.aprox <- 1 - pnorm(q = z))

#> [1] 0.001091522

Note que el \(p_{valor \ aprox.} = 0,0011\) obtenida en la salida anterior, de igual manera, conduce a rechazar la hipótesis nula, dado que es menor que el nivel de significancia establecido. En este caso el \(p_{valor \ aprox.} = 0,0011\) se considera una buena aproximación del \(p_{valor}\) dado que \(np = 9 > 5\).

También, la aproximación del \(p_{valor}\) se puede obtener con la función prop.test del paquete base, como se muestra a continuación, con el siguiente código.

prop.test(
  x = 16, n = 18, p = 0.5, alternative = "greater",
  conf.level = 0.95, correct = TRUE
)

#> 
#>  1-sample proportions test with continuity correction
#> 
#> data:  16 out of 18, null probability 0.5
#> X-squared = 9.3889, df = 1, p-value = 0.001092
#> alternative hypothesis: true p is greater than 0.5
#> 95 percent confidence interval:
#>  0.6803061 1.0000000
#> sample estimates:
#>         p 
#> 0.8888889

Note que la aproximación del \(p_{valor}\) obtenida en estas dos últimas salidas son iguales.

1.1.3 Problemas Propuestos

1. En una población de personas con edades entre 15 y 19 años, se encontró que \(6,9\%\) son positivas a una determinada enfermedad. Supongamos que una muestra aleatoria de 25 de ellas son seleccionadas y que 4 resultaron positivas.
   1. ¿Poseen los datos suficiente evidencias para indicar que la proporción de positivas en esta población es mayor de \(7\%\)? (Use \(\alpha = 5\%\)).
   2. ¿Qué pasa si en vez de tomar la muestra aleatoria de 25, se toma una muestra aleatoria de 30? (Use \(\alpha = 5\%\)).
2. Un fabricante afirma que el \(5\%\) de sus computadores pueden estar defectuosos. Un comprador recibe 20 computadores de los cuales 3 están defectuosos.
   1. ¿Poseen los datos suficiente evidencias para indicar que la proporción de computadoras defectuosas en esta población es menor de \(5\%\)? (Use \(\alpha = 5\%\)).
   2. ¿Qué pasa si en vez de tomar la muestra aleatoria de 20, se toma una muestra aleatoria de 30? (Use \(\alpha = 5\%\)).

1.2 Prueba del Signo

La prueba del signo es un test no paramétrico que se usa para contrastar hipótesis referentes a la mediana de una población. La aplicación de este test requiere que la variable analizada sea cuantitativa, medida al menos en una escala ordinal.

1.2.1 Racionalización de la Prueba del Signo

Supóngase que se tiene una variable aleatoria continua \(X\) con función de distribución de probabilidad acumulada o simplemente distribución de probabilidad F(x) desconocida. Y que se quiere contrastar la hipótesis nula de que la mediana de \(X\) es \(\widetilde{\mu}_{0}\), es decir, \(H_{0}:\widetilde{\mu}= \widetilde{\mu}_{0}\). Si la hipótesis nula es verdadera, la probabilidad de que cualquier valor de \(X\) sea menor o igual que la mediana \(\widetilde{\mu}_{0}\) es \(\frac{1}{2}\), esto es, \(P \left(X< \widetilde{\mu}_{0}\right)=\frac{1}{2}\). Por complemento, la probabilidad de que cualquier valor de \(X\) sea mayor o igual que la mediana también es \(\frac{1}{2}\), es decir, \(P \left(X\geq \widetilde{\mu}_{0}\right)=1-P \left(X\leq \widetilde{\mu}_{0}\right)=1-F \left(\widetilde{\mu}_{0} \right)=\frac{1}{2}\). Sea \(X_1,X_2,\dotsc,X_n\) una muestra de tamaño \(n\) de \(X\), tal que \(X_i\) representa el \(i-ésimo\) valor de \(X\) observado en la muestra, para \(i=1,2,\dotsc,n\). Está claro que los \(X_i\) son variables aleatorias independientes con la misma distribución de \(X\), es decir \(F(x)=F\left(x_i \right)\). Si se define la siguiente variable aleatoria que toma el valor 1, si el \(i-ésimo\) valor muestral es mayor que la mediana y 0 si es menor, lo que se puede expresar de la siguiente manera:

\[ D_{i}=\begin{cases} 1, \ \textit{ si } X_i> \widetilde{\mu}_0\\ 0, \ \textit{ si } X_i< \widetilde{\mu}_0 \end{cases}. \]

De lo anterior se deduce que la probabilidad de encontrar cualquier valor en la muestra que esté por encima de la mediana \(\widetilde{\mu}_0\), de manera teórica, es:

\[ P\left(D_i=1 \right)=P \left(X> \widetilde{\mu}_{0}\right)=1-P \left(X\leq \widetilde{\mu}_{0}\right)=1-F \left(\widetilde{\mu}_{0} \right)=\frac{1}{2}. \]

Mientras que, por complemento, la probabilidad de encontrar un valor cualquiera en la muestra que esté por debajo de la mediana \(\widetilde{\mu}_{0}\) es:

\[ P\left(D_i=0 \right)=P \left(X< \widetilde{\mu}_{0}\right)=P \left(X\leq \widetilde{\mu}_{0}\right)=F \left(\widetilde{\mu}_{0} \right)=\frac{1}{2}. \]

Note que los \(D_i\) constituyen \(n\) réplicas independientes de un experimento de Bernoulli con probabilidad de éxito \(p=P(D_i=1)=\frac{1}{2}\). Ahora, sea \(S\) la variable aleatoria que representa el número de valores muestrales que se encuentran por encima de la mediana propuesta bajo la hipótesis nula \(\widetilde{\mu}_{0}\). Esto es:

\[ \begin{equation} S=\sum_{i=1}^{n}D_i. \end{equation} \tag{1.12} \]

Dado que la probabilidad de que cualquier valor en la muestra se encuentre por encima de la mediana \(\widetilde{\mu}_{0}\) es \(\frac{1}{2}\), y que los \(D_i\) representan \(n\) réplicas independientes de un experimento de Bernoulli, entonces \(S\) es una variable aleatoria binomial con parámetros \(n\) y \(p=\frac{1}{2}\). Por lo tanto su media es \(np=\frac{n}{2}\) y su varianza es \(npq=\frac{n}{4}\). Teóricamente hablando, dado que \(X\) es una variable aleatoria continua, es imposible encontrar en la muestra un \(X_i=\widetilde{\mu}_{0}\), debido a que \(P \left(X=\widetilde{\mu}_{0} \right)= 0\), por ser \(X\) una variable continua. Pero en la práctica esto frecuentemente ocurre, debido generalmente a imprecisión en el instrumento que se ha usado para medir \(X\). Por tal motivo, cuando en la muestra ocurren valores muestrales iguales a la mediana \(\widetilde{\mu}_{0}\), estos deben ser eliminados y el tamaño de la muestra se redefine como el número de observaciones en la muestra para las cuales \(X_i \neq \widetilde{\mu}_{0}\). Como se dijo anteriormente, bajo la hipótesis nula, \(S\sim B\left ( n, p=\frac{1}{2} \right )\), por lo tanto la probabilidad de encontrar de manera exacta \(s\) observaciones en una muestra en particular por encima de \(\widetilde{\mu}_{0}\) (éxitos), viene dada, haciendo \(p=\frac{1}{2}\) en la ecuación (1.2), por:

\[ \begin{equation} P\left ( S=s \right )=p _{B} \left(s, \ n, \ p = \frac{1}{2} \right)= \frac{1}{2^{n}}\binom{n}{s}, \ s = 0, 1, \dotsc, n, \ n>0 \ y \ p =\frac{1}{2}. \tag{1.13} \end{equation} \]

Mientras que la probabilidad de obtener a lo más \(s\) observaciones por encima de \(\widetilde{\mu}_{0}\), se obtiene sustituyendo en la ecuación (1.3) \(p\) por \(\frac{1}{2}\), es decir:

\[ \begin{equation} P(S\leq s)=F_{B}\left (s, \ n, \ p=\frac{1}{2} \right )=\frac{1}{2^{n}}\sum_{i=0}^{s}\binom{n}{i}, \ s=0, 1, \dotsc, n, \ n > 0 \ y \ p=\frac{1}{2}. \tag{1.14} \end{equation} \]

Fíjese que probar la hipótesis nula \(H_{0}:\widetilde{\mu}= \widetilde{\mu}_{0}\) es equivalente a probar la hipótesis nula de que la proporción de posibles valores de \(X\) por encima de \(\widetilde{\mu}_{0}\) es igual a \(\frac{1}{2}\), es decir, \(H_{0}:p= \frac{1}{2}\). En tal sentido, la prueba del signo puede considerarse como un caso particular de la prueba binomial cuando \(p_0=\frac{1}{2}\). Para contrastar la hipótesis nula referente a la proporción se usa el estadístico de prueba S, el cual, bajo la hipótesis nula \(H_{0}:\widetilde{\mu}= \widetilde{\mu}_{0}\), se distribuye binomialmente con parámetros \(n\) y \(p=\frac{1}{2}\), por lo que su media es \(\frac{n}{2}\) y su varianza es \(\frac{n}{4}\). Ahora, si \(s\) es un valor de \(S\) obtenido a partir de una muestra de \(X\), un valor de \(s\) muy alejado de \(\frac{n}{2}\) (por arriba o por debajo), sería evidencia de que \(p \neq \frac{1}{2}\) y en consecuencia \(\widetilde{\mu} \neq \widetilde{\mu}_{0}\).

En este orden de ideas, suponga que se quiere probar la hipótesis alternativa unilateral superior \(H_{1}:\widetilde{\mu}> \widetilde{\mu}_{0}\), lo que sería equivalente a probar la hipótesis alternativa para la proporción \(H_1: p> \frac{1}{2}\), que de acuerdo a la ecuación (1.4), el p.valor sería:

\[ \begin{equation} p_{valor}=1-F_{B}\left ( s-1, \ n, \ p=\frac{1}{2} \right )=1- \frac{1}{2^{n}} \sum_{i=0}^{s-1}\binom{n}{i}. \tag{1.15} \end{equation} \]

En caso de que la hipótesis alternativa para la mediana fuera unilateral inferior \(H_{1}:\widetilde{\mu}< \widetilde{\mu}_{0}\), su hipótesis equivalente para la proporción sería \(H_1: p < \frac{1}{2}\). En tal sentido, de la ecuación (1.5), el \(p_{valor}\) para este caso sería,

\[ \begin{equation} p_{valor}=F_{B}\left (s, \ n, \ p=\frac{1}{2} \right )= \frac{1}{2^{n}} \sum_{i=0}^{s}\binom{n}{i}. \tag{1.16} \end{equation} \]

Para los casos en los que la hipótesis alternativa es unilateral el criterio de decisión consiste en aceptar \(H_{1}\) si el \(p_{valor}< \alpha\) se acepta H_1, ya sea que esta haya sido planteado en forma unilateral superior \(\left(H_{1}:\widetilde{\mu}> \widetilde{\mu}_{0}\right)\) o unilateral inferior \(\left(H_{1}:\widetilde{\mu}< \widetilde{\mu}_{0}\right)\).

Para una hipótesis alternativa de dos colas \(H_{1}:\widetilde{\mu}\neq \widetilde{\mu}_{0}\), su correspondiente en términos de proporción sería \(H_1: p \neq \frac{1}{2}\) , por lo que de la ecuación (1.6), su \(p_{valor}\) sería:

\[ \begin{equation} p_{valor}= \begin{cases} 2\left [ 1-F_{B}\left ( s-1,\ n,\ p = \frac{1}{2} \right ) \right ], \ \textit{si} \ s > \frac{n}{2} \\ 2F_{B}\left ( s,\ n, \ p = \frac{1}{2} \right ), \ \textit{si} \ s \leq \frac{n}{2} \end{cases}. \tag{1.17} \end{equation} \]

R proporciona la función SIGN.test(x, y = NULL, md = 0, alternative = "two.sided", conf.level = 0.95, ...) disponible en el paquete BSDA que permite ejecutar la prueba del signo para una muestra (desarrollada acá) y dos muestras relacionadas (tratada en el siguiente núcleo temático). Los parámetros de esta función se definen en la tabla 1.5.

knitr::kable(
  data.frame(
    stringsAsFactors = FALSE,
    Parámetros = c("x", "y", "md", "alternative", "conf.level", "$\\dotsc$"),
    Descripción = c(
      "Vector numérico; los NA e Inf están permitidos pero serán eliminados.",
      "Vector numérico opcional; los NA e Inf están permitidos pero serán eliminados.",
      "Un número único que represente el valor de la mediana de la población especificada por la hipótesis nula.",
      "Es una cadena de caracteres, que puede tomar el valor: \"greater\", \"less\" o \"two.sided\", o la letra inicial de cada una de ellas, que indica la hipótesis alternativa especificada. Para pruebas de una muestra, alternative se refiere a la verdadera mediana de la población originaria en relación con el valor hipotético de la mediana.",
      "Nivel de confianza para el intervalo de confianza establecido, limitado a estar entre cero y uno.",
      "Otros argumentos que se deben proporcionar a o desde los métodos."
    )
  ),
  booktabs = TRUE,
  caption = "Descripción de los parámetros de la función SIGN.test de R",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 1.5: Descripción de los parámetros de la función SIGN.test de R

Parámetros	Descripción
x	Vector numérico; los NA e Inf están permitidos pero serán eliminados.
y	Vector numérico opcional; los NA e Inf están permitidos pero serán eliminados.
md	Un número único que represente el valor de la mediana de la población especificada por la hipótesis nula.
alternative	Es una cadena de caracteres, que puede tomar el valor: “greater”, “less” o “two.sided”, o la letra inicial de cada una de ellas, que indica la hipótesis alternativa especificada. Para pruebas de una muestra, alternative se refiere a la verdadera mediana de la población originaria en relación con el valor hipotético de la mediana.
conf.level	Nivel de confianza para el intervalo de confianza establecido, limitado a estar entre cero y uno.
\(\dotsc\)	Otros argumentos que se deben proporcionar a o desde los métodos.

De manera análoga a la prueba binomial, cuando el tamaño de la muestra es grande la distribución del estadístico de prueba \(S\) tiende a la normal conforme \(n\) tiende a infinito, en este caso la convergencia es rápida dado que para \(p=\frac{1}{2}\) la distribución binomial es simétrica, por tanto, valores conservadores de \(n\) garantizan una buena aproximación. Luego, dado que bajo la hipótesis nula \(p_{0}=\frac{1}{2}\) y haciendo la sustitución respectiva en la ecuación (1.7), el estadístico \(Z\), queda definido como:

\[ \begin{equation} Z=\frac{2S-n}{\sqrt{n}}. \tag{1.18} \end{equation} \]

Dado que \(Z\) es una estandarización de \(S\), entonces \(Z\simeq N\left(\mu_{Z}=0, \sigma_{Z}=1\right)\) para un \(n\) suficientemente grande. Ahora, tal como se aplicó corrección por continuidad de Yate para conseguir mejores aproximaciones en el cálculo de probabilidades en la prueba binomial dada por la ecuación (1.8), del mismo modo, mejores aproximaciones del \(p_{valor}\) en la prueba del signo se obtienen con el estadístico:

\[ \begin{equation} Z=\frac{\left(2S \pm 1 \right)-n}{\sqrt{n}}, \tag{1.19} \end{equation} \]

Donde \(2S + 1\) se usa cuando \(S < \frac{n}{2}\) y \(2S - 1\) se usa cuando\(S > \frac{n}{2}\).

En este orden de ideas, un valor aproximado del \(p_{valor}\) cuando la hipótesis alternativa es unilateral superior \(\left(H_{1}:\widetilde{\mu}> \widetilde{\mu}_{0}\right)\)), se obtiene sustituyendo en la ecuación (1.9) \(p_{0}\) por \(\frac{1}{2}\). El resultado se presenta en la siguiente ecuación:

\[ \begin{equation} p_{valor}\cong 1-F_{N}\left ( Z= \frac{2s-1-n}{\sqrt{n}}, \ \mu_{Z}=0, \ \sigma_{Z}=1 \right ). \tag{1.20} \end{equation} \]

De manera análoga, cuando la hipótesis alternativa es unilateral inferior \(\left(H_{1}:\widetilde{\mu}< \widetilde{\mu}_{0}\right)\), el \(p_{valor}\) aproximado se obtiene reemplazando en la ecuación (1.10) \(p_{0}\) por \(\frac{1}{2}\), como sigue:

\[ \begin{equation} p_{valor}\cong F_{N}\left ( Z= \frac{2s+1-n}{\sqrt{n}}, \ \mu_{Z}=0, \ \sigma_{Z}=1 \right ). \tag{1.21} \end{equation} \]

Del mismo modo, cuando la prueba de hipótesis es a dos colas o bilateral \(\left(H_{1}:\widetilde{\mu}\neq \widetilde{\mu}_{0}\right)\), el \(p_{valor}\) (aproximado) se obtiene sustituyendo en la ecuación (1.11) \(p_{0}\) por \(\frac{1}{2}\), como sigue:

\[ \begin{equation} p_{valor}\cong \begin{cases} 2\left [ 1-F_{N}\left ( Z=\frac{2s-1-n}{\sqrt{n}}, \ \mu_{Z}=0, \ \sigma_{Z}=1\right ) \right ] , \ \textit{si} \ s> \frac{n}{2} \\ 2 F_{N}\left ( Z= \frac{2s+1-n}{\sqrt{n}} ,\ \mu_{Z}=0, \ \sigma_{Z}=1\right ), \ \textit{si} \ s \leq \frac{n}{2} \end{cases}. \tag{1.22} \end{equation} \]

1.2.2 Implementación en R de la Prueba del Signo

Ejemplo 1.2 (Prueba del Signo) Los siguientes datos constituyen una muestra aleatoria de 15 mediciones de la clasificación de octano de cierto tipo de gasolina:

\[ \begin{array}{cccccccc} \hline 99,0 & 102,3 & 99,8 & 100,5 & 99,7 & 96,2 & 99,1 & 102,5 \\ 103,3 & 97,4 & 100,4 & 98,9 & 98,3 & 98,0 & 101,6 & \\ \hline \end{array} \]

Pruebe la hipótesis nula \(\widetilde{\mu}=98,0\) contra la hipótesis alternativa \(\widetilde{\mu}> 98,0\) en el nivel de significancia 0,01. Este ejemplo es tomado de Johnson (2012, pág. 460).

Para contrastar estas hipótesis primero es conveniente ordenar la muestra de mayor a menor y luego asignar a cada observación un signo positivo (“+”), igual (“=”) o negativo (“-”) según sea mayor, igual o menor que la mediana, respectivamente. Tal como se muestra en la tabla 1.6.

co <- c(
  99.0, 102.3, 99.8, 100.5, 99.7, 96.2, 99.1, 102.5,
  103.3, 97.4, 100.4, 98.9, 98.3, 98.0, 101.6
)
mediana <- 98
alfa <- 0.01
datos2 <- tibble(
  observacion = 1:length(co),
  co = co
) %>%
  mutate(
    signo = as.factor(
      case_when(
        co > mediana ~ "\\+",
        co == mediana ~ "=",
        co < mediana ~ "\\-"
      )
    )
  ) %>%
  arrange(co)
kableExtra::kbl(
  datos2,
  col.names = c("Observación", "Clasificación de octano", "Signo"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  booktabs = TRUE,
  caption = "\\label{tab2:datos2}Clasificación de octano de cierto tipo de gasolina",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  row_spec(which(datos2$signo == "="), color = "red", strikeout = T) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 1.6: Clasificación de octano de cierto tipo de gasolina

Observación	Clasificación de octano	Signo
6	96,2	-
10	97,4	-
14	98,0	=
13	98,3	+
12	98,9	+
1	99,0	+
7	99,1	+
5	99,7	+
3	99,8	+
11	100,4	+
4	100,5	+
15	101,6	+
2	102,3	+
8	102,5	+
9	103,3	+

Note que en la tabla 1.6 la observación 14 es igual que la mediana establecida en la hipótesis nula, por lo que esta observación se elimina de la muestra y \(n\) se establece como 14. La tabla anterior se resume en la tabla 1.7, ejecutando el siguiente código.

conteo_signo <- datos2 %>%
  group_by(signo) %>%
  summarise(conteo = n())
kable(
  x = conteo_signo,
  col.names = c("Signo", "Conteo"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  booktabs = TRUE,
  caption = "\\label{tab2:datos2}Conteo de signos",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2()

Tabla 1.7: Conteo de signos

Signo	Conteo
-	2
+	12
=	1

La hipótesis alternativa que se pretende contrastar es unilateral superior \(H_{1}:\widetilde{\mu}> 98,0\). Para este caso el valor del estadístico de prueba \(s\) y el \(p_{valor}\) vienen dados por las ecuaciones (1.12) y (1.15). El siguiente código permite determinar el valor del estadístico de prueba \(s\) y el \(p_{valor}\), necesarios para contrastar la hipótesis planteada.

s <- as.numeric(conteo_signo %>%
  dplyr::filter(signo == "\\+") %>%
  dplyr::select(conteo))
n <- as.numeric(conteo_signo %>%
  dplyr::filter(signo != "=") %>%
  dplyr::select(conteo) %>%
  sum())
p.valor <- 1 - pbinom(q = s - 1, size = n, prob = 0.5)
paste("El estadístico de prueba  s =", s, "y el p.valor =", p.valor)

#> [1] "El estadístico de prueba  s = 12 y el p.valor = 0.0064697265625"

En la salida anterior se observa que el valor del estadístico de prueba \(s=12\) , mientras que el \(p_{valor} = 0,0065\) . De donde se concluye que existe suficiente evidencia estadística para asegurar que la mediana de la clasificación de octano del combustible es mayor que 98, debido a que el \(p_{valor}\) es menor que el nivel de significancia especificado (\(\alpha =\) 0,01).

Una forma más expedita de ejecutar la prueba del signo en R es a través de la función SING.test del paquete BSDA cuyos parámetros se describen en la tabla 1.5. El siguiente código muestra la ejecución de la prueba del signo por medio de esta función. Note que el \(p_{valor}\) que se muestra en esta salida es igual al obtenido en la salida anterior.

BSDA::SIGN.test(x = co, md = 98, alternative = "g", conf.level = 0.99)

#> 
#>  One-sample Sign-Test
#> 
#> data:  co
#> s = 12, p-value = 0.00647
#> alternative hypothesis: true median is greater than 98
#> 99 percent confidence interval:
#>  98.13627      Inf
#> sample estimates:
#> median of x 
#>        99.7 
#> 
#> Achieved and Interpolated Confidence Intervals: 
#> 
#>                   Conf.Level  L.E.pt U.E.pt
#> Lower Achieved CI     0.9824 98.3000    Inf
#> Interpolated CI       0.9900 98.1363    Inf
#> Upper Achieved CI     0.9963 98.0000    Inf

Dado que la prueba del signo se puede considerar como un caso particular de la prueba binomial cuando \(p_{0}=\frac{1}{2}\), la función binom.test también puede ser usada para realizar la prueba del signo. El siguiente código muestra la ejecución de la prueba del signo a través de esta función. Note que el \(p_{valor}\) obtenido aquí es idénticos a los anteriores.

s <- as.numeric(conteo_signo %>%
  dplyr::filter(signo == "\\+") %>%
  dplyr::select(conteo)) # conflicto con MASS::select
n <- as.numeric(conteo_signo %>% 
  dplyr::filter(signo != "=") %>%
  dplyr::select(conteo) %>% sum())
binom.test(
  x = s, n = n, p = 0.5, alternative = "greater",
  conf.level = 0.99
)

#> 
#>  Exact binomial test
#> 
#> data:  s and n
#> number of successes = 12, number of trials = 14, p-value = 0.00647
#> alternative hypothesis: true probability of success is greater than 0.5
#> 99 percent confidence interval:
#>  0.5217357 1.0000000
#> sample estimates:
#> probability of success 
#>              0.8571429

Como se dijo anteriormente, el \(p_{valor}\) para la prueba del signo, cuando \(n\) es suficientemente grande, se puede aproximar con la ecuación (1.20). El siguiente código permite evaluar tal ecuación. Como se puede ver, el \(p_{valor \ aprox.}\) generado por este procedimiento \(\left(0,0081\right)\) no difiere mucho del \(p_{valor}\) exacto \(\left(0,0065\right)\) obtenido por los métodos anteriores. Note que el \(p_{valor \ aprox.}\) sigue siendo menor que el nivel de significancia establecido para esta prueba, lo que conduce, de igual manera, a rechazar la hipótesis nula \(H_{0}: \widetilde{\mu}=98,0\).

s <- as.numeric(conteo_signo %>%
  dplyr::filter(signo == "\\+") %>%
  dplyr::select(conteo))
n <- as.numeric(conteo_signo %>%
  dplyr::filter(signo != "=") %>%
  dplyr::select(conteo) %>%
  sum())
p.valor.aprox <- 1 - pnorm(q = (2 * s - 1 - n) / sqrt(n))
paste(
  "El p.valor.aprox = ",
  format(
    x = p.valor.aprox,
    decimal.mark = getOption("OutDec")
  )
)

#> [1] "El p.valor.aprox =  0.008078466"

Otra manera de calcular el \(p_{valor \ aprox.}\) es usando la función prop.test, de manera análoga a la prueba binomial.

s <- as.numeric(conteo_signo %>%
  dplyr::filter(signo == "\\+") %>%
  dplyr::select(conteo))
n <- as.numeric(conteo_signo %>%
  dplyr::filter(signo != "=") %>%
  dplyr::select(conteo) %>%
  sum())
prop.test(
  x = s, n = n, p = 0.5, alternative = "greater",
  conf.level = 0.99, correct = TRUE
)

#> 
#>  1-sample proportions test with continuity correction
#> 
#> data:  s out of n, null probability 0.5
#> X-squared = 5.7857, df = 1, p-value = 0.008078
#> alternative hypothesis: true p is greater than 0.5
#> 99 percent confidence interval:
#>  0.5106275 1.0000000
#> sample estimates:
#>         p 
#> 0.8571429

1.2.3 Problemas Propuestos

1. Se toman 10 muestras de un baño de cultivo sobre placa utilizado en un proceso de fabricación de componentes electrónicos, y se mide el pH del baño. Los valores de pH medidos son 7,91; 7,85; 6,82; 8,01; 7,46; 6,95; 7,05; 7,35; 7,25; 7,42. Los ingenieros plantean que el pH neutro es 7. ¿La muestra indica que esta proposición es correcta? Use \(\alpha = 5\%\).

2. Los siguientes datos representan el número de horas que un temporizador opera antes de que deba recargarse: 1,5; 2,2; 0,9; 1,3; 2,0; 1,6; 1,8; 1,5; 2,0; 1,2 y 1,7. Utiliza la prueba de los signos para probar la hipótesis al nivel de significancia de 0,05 que este temporizador en particular opera con una mediana de 1,8 horas antes de requerir una recarga.

1.3 Prueba de Rangos Asignados de Wilcoxon

La prueba de rangos asignados de Wilcoxon, al igual que la prueba del signo discutida en la sección 1.2, se usa para contrastar la hipótesis de que la mediana de una población es igual a un valor especificado, es decir, \(H_{0}:\widetilde{\mu}= \widetilde{\mu}_{0}\). Esta prueba es más potente que la prueba del signo dado que considera, además del signo de las diferencias entre cada valor observado en la muestra y la mediana propuesta en la hipótesis nula \(X_i-\widetilde{\mu}_{0}\), la magnitud de esta diferencia. Por lo que este test es considerado como la opción más apropiada a la prueba de la mediana cuando no se cumple el supuesto de normalidad de la población bajo estudio.

1.3.1 Racionalización de la Prueba de Rangos Asiganados de Wilcoxon

Suponga que se tiene una variable aleatoria \(X\) cuantitativa medida al menos en una escala ordinal y que \(X_1, X_2, \dotsc, X_n\) es cualquier muestra aleatoria de tamaño \(n\) de \(X\). Wilcoxon (1945) propuso el estadístico

\[ \begin{equation} T^+=\sum_{i=1}^{n}S_iR_i. \tag{1.23} \end{equation} \]

El cual se conoce como estadístico de rangos asignados de Wilcoxon. Donde,

\[ \begin{equation} S_i= \begin{cases} 1, \quad si \; D_i>0\\ 0, \quad si \; D_i<0 \end{cases}. \tag{1.24} \end{equation} \]

Siendo \(D_i= X_i-\widetilde{\mu}_{0}\) y \(R_i=r \left(|D_i| \right)\). En donde la función \(r(.)\) determina los ranqueos, en este caso, ranquea los valores absolutos de los \(D_i\). Esta función se ejecuta con la función rank de la distribución base de R.

Wilcoxon determinó la distribución exacta del estadístico \(T^+\) asumiendo que si \(x_1, x_2, \dotsc,x_n\) representa una muestra de una variable aleatoria continua \(X\), \(x_i \neq \widetilde{\mu}_{0}\), para \(i=1, 2, \dotsc, n\), es decir los valores muestrales eran diferente a la mediana propuesta en la hipótesis nula. Además, si \(\left(|d_r| \right)\) y \(\left(|d_s| \right)\) son los valores absolutos de las diferencias de las observaciones \(r\) y \(s\) a la mediana propuesta en la hipótesis nula \(\widetilde{\mu}_{0}\). Esto es, \(\left(|d_r| \right) \neq \left(|d_s| \right)\), para \(r \neq s\).

Las funciones:

dsignrank(x, n, log = FALSE)
psignrank(q, n, lower.tail = TRUE, log.p = FALSE)
qsignrank(p, n, lower.tail = TRUE, log.p = FALSE)
rsignrank(nn, n)

del paquete stats permiten evaluar la función de densidad, la función de distribución, la función de cuantiles y generar números aleatorios, respectivamente, de la distribución de rangos asignados de Wilcoxon.

Bajo las suposiciones antes descritas, Wilcoxon demostró que la media de \(T^+\) es,

\[ \begin{equation} E \left(T^+| \widetilde{\mu}= \widetilde{\mu}_{0} \right)=\frac{n(n+1)}{4}; \tag{1.25} \end{equation} \]

y la varianza está determinada por

\[ \begin{equation} V \left(T^+| \widetilde{\mu}= \widetilde{\mu}_{0} \right)=\frac{n(n+1)(2n+1)}{24}. \tag{1.26} \end{equation} \]

Suponga que se quiere contrastar la hipótesis nula \(H_0 = \widetilde{\mu}= \widetilde{\mu}_{0}\) contra la hipótesis alternativa \(H_0 = \widetilde{\mu}> \widetilde{\mu}_{0}\). Entonces el \(p_{valor}\) viene dado por,

\[ \begin{equation} p_{valor}=P \left(T^+ \ge t^+ \right)=1 - P \left(T^+ \le t^+ -1 \right). \end{equation} \]

Como \(T^+\) tiene distribución de rangos asignados de Wilcoxon con parámetro \(n\), lo cual se denota como \(T^+ \sim W \left( n \right)\), la ecuación anterior puede reescribirse como:

\[ \begin{equation} p_{valor} = 1 - F_W \left(t^+ - 1, \; n \right). \tag{1.27} \end{equation} \]

Si la hipótesis que se quiere contrastar es que \(H_0: \widetilde{\mu}= \widetilde{\mu}_{0}\) versus la hipótesis alternativa \(H_0: \widetilde{\mu} < \widetilde{\mu}_{0}\). El \(p_{valor}\) viene dado por,

\[ \begin{equation} p_{valor}=P \left(T^+ \le t^+ \right)=P \left(T^+ \le t^+ \right). \end{equation} \]

En resumen, el \(p_{valor}\) cuando la prueba de hipótesis es unilateral inferior, viene dada por:

\[ \begin{equation} p_{valor} = F_W \left(t^+ , \; n \right). \tag{1.28} \end{equation} \]

Ahora, si la prueba de hipótesis es bilateral. Es decir, \(H_0: \widetilde{\mu}= \widetilde{\mu}_{0}\) \(H_0: \widetilde{\mu} \neq \widetilde{\mu}_{0}\). El \(p_{valor}\) (aproximado) viene dado por:

\[ \begin{equation} p_{valor}= \begin{cases} 2\left [ 1-F_{W}\left ( t^+ - 1,\; n \right ) \right ], \; \textit{si} \; t^+ > \frac{n \left(n+1 \right)}{4} \\ 2F_{W}\left ( t^+,\; n \right ), \; \textit{si} \; t^+ \leq \frac{n \left(n+1 \right)}{4} \end{cases}. \tag{1.29} \end{equation} \]

Por otro lado, dado que \(T^+\) es una suma de variables aleatorias, por el teorema del límite central, para un tamaño de muestra grande (\(n\geq 15\)), el estadístico:

\[ \begin{equation} Z= \frac{T^+ \pm 0,5-\frac{n(n+1)}{4}}{\sqrt{\frac{n(n+1)(2n+1)}{24}}} \simeq N \left(\mu=0, \sigma=1\right). \tag{1.30} \end{equation} \]

En los casos donde \(x_i=\widetilde{\mu}_{0}\) se procede a eliminar estas observaciones de la muestra y se toma \(n\) como el número de observaciones en las que \(x_i \neq \widetilde{\mu}_{0}\).

Para una prueba unilateral superior, el \(p_{valor}\), usando el estadístico dado en la ecuación (1.30), se determina con,

\[ \begin{equation} p_{valor} = 1 - F_N \left(Z= \frac{T^+ - 0,5-\frac{n(n+1)}{4}}{\sqrt{\frac{n(n+1)(2n+1)}{24}}}; \; \mu = 0, \; \sigma = 1 \right). \tag{1.31} \end{equation} \]

Si la prueba es unilateral inferior, cuando se usa el estadístico dado en la ecuación (1.30), el \(p_{valor}\) se aproxima con:

\[ \begin{equation} p_{valor} = F_N \left(Z= \frac{T^+ + 0,5-\frac{n(n+1)}{4}}{\sqrt{\frac{n(n+1)(2n+1)}{24}}}; \; \mu = 0, \; \sigma = 1 \right). \tag{1.32} \end{equation} \]

Ahora, si la prueba es a dos colas, el \(p_{valor}\) se calcula con:

\[ \begin{equation} p_{valor}= \begin{cases} 2\left [ 1-F_{N}\left ( \frac{t^+ - 0,5 - \frac{n(n + 1)}{4}}{\sqrt{\frac{n(n+)(2n +1)}{24}}}; \; \mu_Z = 0, \sigma_Z = 1 \right ) \right ], \; \textit{si} \; t^+ > \frac{n \left(n+1 \right)}{4} \\ F_{N}\left ( \frac{t^+ + 0,5 - \frac{n(n + 1)}{4}}{\sqrt{\frac{n(n+)(2n + 1)}{24}}}; \; \mu_Z = 0, \sigma_Z = 1 \right ), \; \textit{si} \; t^+ \le \frac{n \left(n+1 \right)}{4} \end{cases}. \tag{1.33} \end{equation} \]

Cuando la muestra es grande(\(n \ge 15\)) y ocurren diferencias iguales en valor absoluto, es decir, \(\left(|d_r| \right) = \left(|d_s| \right)\), con \(r \ne s\), se asigna el promedio de los ranqueos a los grupos empatados y se realiza una corrección a la varianza (1.26). Esta corrección fue propuesta por Pratt (1.959), como se muestra a continuación:

\[ \begin{equation} V \left(T^+| \widetilde{\mu}= \widetilde{\mu}_{0} \right)=\frac{n(n+1)(2n+1)}{24}-\sum_{j=1}^{g} \frac{t_i \left(t_i^2-1 \right)}{48}. \tag{1.34} \end{equation} \]

Donde \(g\) es el número de grupos y \(t_j\) es el número de observaciones empatadas en el grupo \(j\). Como se puede notar, una observación no empatada es considerada como un grupo de tamaño 1. En particular, si no hay valores de \(\left|D_i \right|\) iguales, entonces \(g = n\) y \(t_j = 1\) para \(j = 1, 2, \dotsc, n\). En este caso, cada término en la ecuación (1.34) de la forma \(t_j (t_j^2 − 1)\) se reduce a cero, y la varianza dada en la ecuación (1.34) se reduce a la varianza del estadístico \(T^+\) cuando no hay empates, mostrada en la ecuación (1.26).

De lo anterior, el estadístico dado en la ecuación (1.30) se redefine de la siguiente manera,

\[ \begin{equation} Z= \frac{T^+ \pm 0,5-\frac{n(n+1)}{4}}{\sqrt{\frac{n(n+1)(2n+1)}{24}-\sum_{j=1}^{g} \frac{t_i \left(t_i^2-1 \right)}{48}}} \simeq N \left(\mu=0, \sigma=1\right). \tag{1.35} \end{equation} \]

Cuando se usa el estadístico anterior para pruebas de hipótesis unilateral superior, el \(p_{valor}\) se obtiene evaluando la siguiente expresión:

\[ \begin{equation} p_{valor} = 1 - F_N \left(Z= \frac{T^+ - 0,5-\frac{n(n+1)}{4}}{\sqrt{\frac{n(n+1)(2n+1)}{24}-\sum_{j=1}^{g} \frac{t_i \left(t_i^2-1 \right)}{48}}}; \; \mu_Z = 0, \; \sigma_Z = 1 \right). \tag{1.36} \end{equation} \]

Para una prueba unilateral inferior usando el estadístico dado en la ecuación (1.35), el \(p_{valor}\) se obtiene con:

\[ \begin{equation} p_{valor} = F_N \left(Z= \frac{T^+ + 0,5-\frac{n(n+1)}{4}}{\sqrt{\frac{n(n+1)(2n+1)}{24}-\sum_{j=1}^{g} \frac{t_i \left(t_i^2-1 \right)}{48}}}; \; \mu_Z = 0, \; \sigma_Z = 1 \right). \tag{1.36} \end{equation} \]

Para un prueba bilateral, usando el estadístico mostrado en la ecuación (1.35), el \(p_{valor}\) se halla usando:

\[ \begin{equation} p_{valor}= \begin{cases} 2\left [ 1-F_{N}\left ( \frac{t^+ - 0,5 - \frac{n(n + 1)}{4}}{\sqrt{\frac{n(n+1)(2n+1)}{24}-\sum_{j=1}^{g} \frac{t_i \left(t_i^2-1 \right)}{48}}}; \; \mu_Z = 0, \sigma_Z = 1 \right ) \right ], \; \textit{si} \; t^+ > \frac{n \left(n+1 \right)}{4} \\ F_{N}\left ( \frac{t^+ + 0,5 - \frac{n(n + 1)}{4}}{\sqrt{\frac{n(n+1)(2n+1)}{24}-\sum_{j=1}^{g} \frac{t_i \left(t_i^2-1 \right)}{48}}}; \; \mu_Z = 0, \sigma_Z = 1 \right ), \; \textit{si} \; t^+ \le \frac{n \left(n+1 \right)}{4} \end{cases}. \tag{1.33} \end{equation} \]

En Conover (1999, pág. 353) se propone el siguiente estadístico de prueba cuando el tamaño de muestra es grande.

\[ \begin{equation} T =T^+ - T^-=2T^+-\frac{n(n+1)}{2}. \tag{1.37} \end{equation} \]

Donde,

\[ T^-=\sum_{i=1}^{n} \left(1- S_i \right)R_i. \]

Por el teorema del límite central la distribución del estadístico \(T\) converge asintóticamente a la distribución normal conforme \(n \rightarrow \infty\). La media de este estadístico, bajo la hipótesis nula, es 0, es decir, \(E \left(T \;| \; \widetilde{\mu}= \widetilde{\mu}_{0} \right) =0\); mientras que su varianza, suponiendo que no es posible encontrar observaciones empatadas en la muestra es

\[ \begin{equation} V \left(T \;| \; \widetilde{\mu}= \widetilde{\mu}_{0} \right)=\frac{n(n+1)(2n+1)}{6}. \end{equation} \]

De lo anterior se obtiene que,

\[ \begin{equation} Z= \frac{T \pm 0,5}{\sqrt{\frac{n(n+1)(2n+1)}{6}}} \simeq N \left(\mu=0, \sigma=1\right). \tag{1.38} \end{equation} \]

Para el estadístico \(T\) dado en la ecuación (1.37) el \(p_{valor}\) para una prueba unilateral superior cuando no hay empates, viene dado por:

\[ \begin{equation} p_{valor} = 1 - F_N \left(Z= \frac{t - 0,5}{\sqrt{\frac{n(n+1)(2n+1)}{6}}}; \; \mu_Z = 0, \; \sigma_Z = 1 \right). \tag{1.31} \end{equation} \]

Cuando la prueba es unilateral inferior y se usa el estadístico dado en la ecuación (1.37) sin presencia de empates, el \(p_{valor}\) es:

\[ \begin{equation} p_{valor} =F_N \left(Z= \frac{t + 0,5}{\sqrt{\frac{n(n+1)(2n+1)}{6}}}; \; \mu_Z = 0, \; \sigma_Z = 1 \right). \tag{1.31} \end{equation} \]

Ahora, si la prueba es bilateral y no hay empates, el \(p_{valor}\) se calcula con,

\[ \begin{equation} p_{valor}= \begin{cases} 2\left [ 1-F_{N}\left ( \frac{t - 0,5 }{\sqrt{\frac{n(n+)(2n + 1)}{6}}}; \; \mu_Z = 0, \sigma_Z = 1 \right ) \right ], \; \textit{si} \; t > 0 \\ F_{N}\left ( \frac{t + 0,5 }{\sqrt{\frac{n(n+)(2n + 1)}{6}}}; \; \mu_Z = 0, \sigma_Z = 1 \right ), \; t \le 0 \end{cases}. \tag{1.39} \end{equation} \]

Cuando hay empates en la muestra y esta es grande (\(n \ge 15\)), el estadístico dado en la ecuación (1.37) tiene \(E \left(T \;| \; \widetilde{\mu}= \widetilde{\mu}_{0} \right) =0\); y varianza

\[ \begin{equation} V \left(T \;| \; \widetilde{\mu}= \widetilde{\mu}_{0} \right)=\sum_{i=1}^{n}R_i^2. \end{equation} \]

En consecuencia,

\[ \begin{equation} Z= \frac{T \pm 0,5}{\sqrt{\sum_{i=1}^{n}R_i^2}} \simeq N \left(\mu=0, \sigma = 1\right). \tag{1.40} \end{equation} \]

Por lo que, para muestras grandes y sin empates, haciendo uso del estadístico dado en la ecuación anterior, el calculo del \(p_{valor}\) para pruebas de hipótesis unilateral superior, unilateral inferior y bilateral viene dadas por las siguientes ecuaciones, respectivamente.

\[ \begin{equation} p_{valor} = 1 - F_N \left(Z= \frac{t - 0,5}{\sqrt{\sum_{i=1}^{n}R_i^2}}; \; \mu_Z = 0, \; \sigma_Z = 1 \right). \tag{1.41} \end{equation} \]

\[ \begin{equation} p_{valor} =F_N \left(Z= \frac{t + 0,5}{\sqrt{\sum_{i=1}^{n}R_i^2}}; \; \mu_Z = 0, \; \sigma_Z = 1 \right). \tag{1.41} \end{equation} \]

\[ \begin{equation} p_{valor}= \begin{cases} 2\left [ 1-F_{N}\left ( \frac{t - 0,5 }{\sqrt{\sum_{i=1}^{n}R_i^2}}; \; \mu_Z = 0, \sigma_Z = 1 \right ) \right ], \; \textit{si} \; t > 0 \\ 2F_{N}\left ( \frac{t + 0,5 }{\sqrt{\sum_{i=1}^{n}R_i^2}}; \; \mu_Z = 0, \sigma_Z = 1 \right ), \; t \le 0 \end{cases}. \tag{1.42} \end{equation} \]

1.3.2 Implementación en R de la Prueba de Rangos Asignados de Wilcoxon

Ejemplo 1.3 (Prueba de Rangos Asignados de Wilcoxon) Use la prueba de rangos asignados de Wilcoxon a los datos del ejemplo 1.2 para contrastar las hipótesis planteadas en el ejemplo.

La siguiente tabla muestra los rangos asignados para la concentración de octano.

co <- c(
  99.0, 102.3, 99.8, 100.5, 99.7, 96.2, 99.1, 102.5,
  103.3, 97.4, 100.4, 98.9, 98.3, 98.0, 101.6
)
mediana <- 98
tabla_co <- tibble(
  i = 1:length(co),
  co = co,
  dif = co - mediana,
  va_dif = na_if(abs(dif), 0),
  rango = sign(dif) * rank(
    va_dif, ties.method = "average") %>% 
    na_if(max(.))
) %>% 
  arrange(va_dif)
options(knitr.kable.NA = "")
kableExtra::kbl(
  tabla_co,
  col.names = c(
    "$i$", "Concentración de octano", "$D_i$",  
    "$\\left| D_i \\right|$","$R_i$"
  ),
  format.args = list(decimal.mark = ",", big.mark = "."),
  booktabs = TRUE,
  caption = "\\label{tab2:datos-co}Ranqueo de la concentración de octano",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  row_spec(which(tabla_co$dif == 0), color = "red",
           strikeout = TRUE) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 1.8: Ranqueo de la concentración de octano

\(i\)	Concentración de octano	\(D_i\)	\(\left| D_i \right|\)	\(R_i\)
13	98,3	0,3	0,3	1,0
10	97,4	-0,6	0,6	-2,0
12	98,9	0,9	0,9	3,0
1	99,0	1,0	1,0	4,0
7	99,1	1,1	1,1	5,0
5	99,7	1,7	1,7	6,0
3	99,8	1,8	1,8	7,5
6	96,2	-1,8	1,8	-7,5
11	100,4	2,4	2,4	9,0
4	100,5	2,5	2,5	10,0
15	101,6	3,6	3,6	11,0
2	102,3	4,3	4,3	12,0
8	102,5	4,5	4,5	13,0
9	103,3	5,3	5,3	14,0
14	98,0	0,0

En el siguiente script se muestra el calculo del estadístico \(T^+\) y el \(p_{valor}\) para una prueba unilateral superior, dados en las ecuaciones (1.23) y (1.27).

Tmas <- tabla_co %>%
  dplyr::select(
    rango
  ) %>%
  dplyr::filter(
    rango > 0
  ) %>%
  sum(., na.rm = TRUE)
n <- nrow(drop_na(tabla_co))
p.valor <- 1 - psignrank(
  q = Tmas - 1, n = n, lower.tail = TRUE
)
paste("El valor del estadístico Tmas =", Tmas)
paste(
  "El p.valor =", format(
    round(p.valor, digits = 6),
    mode = "double",
    big.mark = ".",
    decimal.mark = ","
  )
)

#> [1] "El valor del estadístico Tmas = 95.5"
#> [1] "El p.valor = 0,002014"

Como se muestra en la salida anterior, el estadístico \(T^+ =95,5\). Mientras que el \(p_{valor} = 0,002014\), lo cual implica que hay suficiente evidencia estadística para asegurar con un 99% de que la mediana de la concentración de octano es superior a 98.

Note que el \(p_valor\) obtenido con la prueba del signo (0,0065) en el ejemplo 1.2 no difiere mucho del obtenido con esta prueba (0,002014). En ambos casos se rechaza la hipótesis nula, dado que \(p_{valor} < 0,01\).

La prueba de rangos con signos de Wilcoxon se puede aplicar, de manera directa, con la función wilcoxon.test del paquete exactRankTests.

wilcox.exact(
  x = co, y = NULL, alternative = "g", mu = 98,
  paired = FALSE, exact = TRUE, conf.int = TRUE,
  conf.level = 0.99
)

#> 
#>  Exact Wilcoxon signed rank test
#> 
#> data:  co
#> V = 95.5, p-value = 0.002075
#> alternative hypothesis: true mu is greater than 98
#> 99 percent confidence interval:
#>  98.35   Inf
#> sample estimates:
#> (pseudo)median 
#>         99.775

Dado que el tamaño de la muestra es razonablemente grande el \(p_{valor}\) puede aproximarse haciendo uso de la convergencia asintótica a la normal del estadístico \(T^+\) con la ecuación (1.31). El cálculo se muestra en el siguiente trozo de código.

media <- n*(n + 1)/4
var <- n*(n + 1)*(2*n + 1)/24
z <- (Tmas - 0.5 - media)/sqrt(var)
p.valor <- 1 - pnorm(q = z)
paste("El valor del estadístico Z =", round(x = z, digits = 4))
paste("El p.valor =", round(x = p.valor, digits = 4))

#> [1] "El valor del estadístico Z = 2.668"
#> [1] "El p.valor = 0.0038"

Como se observa, en la tabla 1.8 hay presencia de ceros y rangos empatados por lo que conviene calcular el \(p_{valor}\) usando la ecuación (1.36). Para ello es necesario construir la tabla de rango empatados, lo cual se muestra con el siguiente script.

empates <- tabla_co %>%
  dplyr::select(rango) %>% 
  drop_na(.) %>%
  dplyr::transmute(
    rango = abs(rango)
  ) %>% 
  table(.) %>% 
  as.data.frame(.) %>%
  dplyr::mutate(
    j = 1:length(Freq)
  ) %>% 
  dplyr::rename(
    rango = ".",
    tj = "Freq"
  ) %>% 
  dplyr::select(c(3, 1, 2))
kableExtra::kbl(
  empates,
  col.names = c(
    "$j$", "Rango", "$T_j$"
  ),
  format.args = list(decimal.mark = ",", big.mark = "."),
  booktabs = TRUE,
  caption = "\\label{tab2:empates-co}Rangos empatados",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 1.9: Rangos empatados

\(j\)	Rango	\(T_j\)
1	1	1
2	2	1
3	3	1
4	4	1
5	5	1
6	6	1
7	7.5	2
8	9	1
9	10	1
10	11	1
11	12	1
12	13	1
13	14	1

El calculo del \(p_{valor}\), cuando hay empates, se obtiene con el siguiente script.

cor_var <- sum(empates$tj * (empates$tj^2 - 1)) / 48
var.cor <- var - cor_var
z.cor <- (Tmas - 0.5 - media)/sqrt(var.cor)
p.valor <- 1 - pnorm(q = z.cor)
paste("El valor del estadístico Z =", round(x = z.cor, digits = 4))
paste("El p.valor =", round(x = p.valor, digits = 6))

#> [1] "El valor del estadístico Z = 2.6687"
#> [1] "El p.valor = 0.003808"

La función wilcox.test del paquete stats también realiza esta prueba, como se muestra en el siguiente código.

wilcox.test(
  x = co, y = NULL, alternative = "g", mu = 98,
  paired = F, exact = T, conf.int = T, conf.level = 0.99
)

#> 
#>  Wilcoxon signed rank test with continuity correction
#> 
#> data:  co
#> V = 95.5, p-value = 0.003808
#> alternative hypothesis: true location is greater than 98
#> 99 percent confidence interval:
#>  98.30008      Inf
#> sample estimates:
#> (pseudo)median 
#>       99.80006

Nota: La función wilcox.test no ejecuta la prueba exacta en presencia de ceros y/o empates.

Otra forma de hallar el \(p_{valor}\), cuando hay empates, es usando la convergencia asintótica a la normal del estadístico \(T\), por medio de la ecuación (1.41).

T <- sum(tabla_co$rango, na.rm = TRUE)
var_T <- sum(tabla_co$rango ^ 2, na.rm = TRUE)
z <- (Tmas - 0.5)/sqrt(var_T)
p.valor <- 1 - pnorm(q = z)
paste("El valor del estadístico Z =", round(x = z, digits = 6))
paste("El p.valor =", round(x = p.valor, digits = 6))

#> [1] "El valor del estadístico Z = 2.982618"
#> [1] "El p.valor = 0.001429"

1.4 Prueba de Bondad de Ajuste Ji-cuadrado

La prueba Ji-cuadrado de bondad de ajuste se usa para contrastar la hipótesis de que una variable aleatoria tiene una determinada distribución de probabilidad. Para ello la variable aleatoria debe ser categórica, discreta o continua (preferiblemente categórica o discreta). En este último caso la variable aleatoria debe ser categorizada, para ello el recorrido de la variable debe ser particionado en \(k\) intervalos (construir tabla de distribución de frecuencia en intervalos de clases).

1.4.1 Racionalización de la Prueba de Bondad de Ajuste Ji-cuadrado

Sea \(X\) una variable aleatoria con distribución de probabilidad \(F\left(x,\theta\right)\), donde \(\theta=\left(\theta_1, \dotsc, \theta_r\right)\) es un vector de dimensión \(r\) que define los parámetros de la distribución; si \(\left(x_1,\dotsc, x_n \right)\) es una m.a.s. de \(X\); entonces la prueba de bondad de ajuste Ji-cuadrado permite contrastar, a la luz de la información de la muestra, si la distribución de \(X\) es la indicada en la hipótesis nula \(\left(H_0: F(x,\theta)=F(x,\theta_0 )\right)\). Se van a distinguir dos casos, según que \(F\) sea perfectamente conocida, es decir, el vector de parámetros bajo la hipótesis nula, \(\theta_0=\left(\theta_{10}, \dotsc, \theta_{r0} \right)\), es conocido por ejemplo viene dada por la normal de media 0 y de desviación típica 1, o bien que \(F\) sea conocida pero se desconozca \(\theta_0\), en tal caso \(\theta_0\) debe ser estimado con la información de la muestra; por ejemplo, que \(F\) venga dada por la normal de media \(\mu\) y de desviación típica \(\sigma\) desconocidas, en esta situación \(\mu\) y \(\sigma\) deberán ser estimados con la información muestral.

En el primer caso, suponga que \(I\) define el recorrido de la población \(X\) y que \(I\) se ha dividido en \(k\) intervalos mutuamente excluyentes y colectivamente exhaustivos, en principio en forma arbitraria, \(I_1, \dotsc, I_k\) con lo que es posible calcular bajo la hipótesis nula las probabilidades \(p_{i0}=P\left\{x\in I_i\right\}\), para \(i = 1, \dotsc, k\).

Si la muestra es compatible con \(F\left(x,\theta_0 \right)\) se puede suponer la hipótesis nula \(H_0:\) consistente en que el vector \((n_1, \dotsc, n_k )\), con \(n_i\) número de veces que en la muestra aparece un elemento del intervalo \(I_i\), tiene una distribución \(\left (N_{1} = n_{1}, \dotsc, N_{k} = n_{k} \right )\sim \textit{Multinomial}\left ( n, p_{10}, \dotsc, p_{k0} \right )\), frente a la hipótesis alternativa \(H_1:\left (N_{1} = n_{1}, \dotsc, N_{k} = n_{k} \right )\sim \textit{Multinomial}\left ( n, p_{1}, \dotsc, p_{k} \right )\) . La única diferencia está en que bajo la hipótesis nula los parámetros de la distribución multinomial \(p_0 =(p_{10},\dotsc,p_{k0})\) son conocidos, mientras que bajo la hipótesis alternativa los parámetros \(p = (p_1,\dotsc,p_k )\) son desconocidos.

Si se emplea el contraste de la razón de verosimilitudes, se tiene que

\[ \lambda\left ( n_{1}, \dotsc, n_{k} \right ) =\frac{\frac{n!}{n_1! \dotsm n_k!}\left ( p_{10} \right )^{n_1} \dotsm \left ( p_{k0} \right )^{n_k}}{\underset{p_1, \dotsc, p_k}{SUP}\frac{n!}{n_1! \dotsm n_k!}\left ( p_1 \right )^{n_1} \dotsm \left ( p_k \right )^{n_k}}. \]

El supremo del denominador, asumiendo que \(\sum_{i=1}^{k}p_i=1\), se alcanza en \(\widehat{p}_i=\frac{n_i}{n}\), dado que estos son los estimadores de máxima verosimilitud de los \(p_i\). Por tanto, la región crítica será de la forma

\[ RC=\left\{\lambda \left ( n_1, \dotsc, n_k \right ) = \prod_{i=1}^{k} \left ( \frac{p_{i0}}{\widehat{p}_i} \right )^{n_i} \leq k \right\}. \]

Bajo \(H_0\) la distribución asintótica de \(-2\ln\lambda(n_1, \dotsc, n_k )\) es una \(\chi^{2}\) con \(k-1\) grados de libertad, ya que, \(H_0∪H_1\) tiene de dimensión \(k-1\) y la hipótesis nula tiene de dimensión 0. Para un tamaño \(\alpha\), se obtiene la región crítica

\[ RC = \left\{-2 \ln \lambda \left ( n_1, \dotsc, n_k \right )=-2\sum_{i=1}^{k}n_i \left [ \ln p_{i0}- \ln\widehat{p}_i\right ]\geq \chi ^{2}_{k-1};1-\alpha \right\}. \]

K. Pearson introdujo un procedimiento alternativo de resolver el problema, mediante una distancia ponderada entre el \(EMV\) del vector multinomial \(\widehat{p}\) y el valor \(p_0\) supuesto bajo la hipótesis nula, en concreto la cantidad

\[ D^{2} \left ( p_0, \widehat{p} \right ) = \sum_{i=1}^{k}\lambda_i\left ( \widehat{p}_i -p_{i0} \right )^{2} \]

con

\[ \lambda_i=\frac{n}{n_i} \]

y rechazar la hipótesis nula para valores grandes de la distancia; por lo que la región crítica es

\[ RC=\left\{ \sum_{i=1}^{k}\frac{\left ( n_i - np_{i0} \right )^{2}}{np_{i0}}\geq c\right\} \]

Con el fin de poder determinar la constante \(c\) para que el contraste tenga un nivel de significancia de tamaño \(\alpha\), es necesario conocer la distribución de la cantidad

\[ \sum_{i=1}^{k}\frac{\left ( n_i - np_{i0} \right )^{2}}{np_{i0}} \]

bajo la hipótesis nula, resultando que asintóticamente tiene la misma distribución que \(-2\ln\lambda(n_1,\dotsc,n_k )\). Una justificación de este hecho radica en observar que

\[ -2\ln\lambda(n_1,\dotsc,n_k )=-2\sum_{i=1}^{k}n_i \ln\left( \frac {p_{i0}}{\widehat{p}_i} \right) \]

y si se emplea que, aproximadamente, \(\ln(x)≃(x-1)-\frac{1}{2} (x-1)^2\), se tiene que

\[ \begin{equation} \begin{split} -2\ln\lambda(n_1,\dotsc,n_k ) & \simeq -2\sum_{i=1}^{k}n_i \left[ \left( \frac{p_{i0}}{\widehat{p}_i}- 1 \right) - \frac{1}{2} \left( \frac{p_{i0}}{\widehat{p}_i} -1 \right)^{2} \right] \\ & \simeq \sum_{i=1}^{k}n_i \left( \frac{p_{i0}}{\widehat{p}_i} -1 \right)^{2} \\ & \simeq \sum_{i=1}^{k} \frac{\left( n_{i} - np_{i0} \right)^2}{n_{i}} \\ & \simeq \sum_{i=1}^{k} \frac{\left( n_{i} -np_{i0} \right)^2}{np_{i0}}, \end{split} \end{equation} \]

donde la última aproximación se debe a que \(n_i\simeq np_{i0}\). Por lo que de manera resumida

\[ \begin{equation} D^2 (p_0,\widehat{p}_i )=\sum_{i=1}^{k} \frac{\left( n_{i} -np_{i0} \right)^2}{np_{i0}} \xrightarrow{cs }\chi_{k-1}^{2}. \tag{1.43} \end{equation} \]

Como se ha mostrado, la distribución del estadístico dado en la ecuación (1.43), bajo la hipótesis nula \(\left(p_i=p_{i0} \right)\), tiende hacia la distribución ji-cuadrado con \(k-1\) grados de libertad cuando \(n\) es grande. Lo grande que debe ser \(n\), para que la aproximación sea aceptablemente buena, está en relación con los valores de las probabilidades \(p_i\); usualmente se considera que \(np_i≥5 \ \textit{para cada} \;i=1,2,\dotsc,k\).

Ello sirve, a menudo, para delimitar los conjuntos \(I_i\), que han de elegirse de forma que se cumplan las restricciones \(np_i≥5\).

Sin embargo, puesto que el contraste no discrimina entre distribuciones que asignen la misma probabilidad a los conjuntos \(I_i\), es aconsejable que el número \(k\) de particiones no sea inferior a 5 (salvo en ajustes a distribuciones discretas o cualitativas con menor número de valores posibles). Así, tendrá que ser alguna de las \(p_{i0}≤ \frac{1}{5}\) y, por consiguiente, es necesario que sea \(n≥25\). Aunque, evidentemente, se conseguirá mayor exactitud cuanto más grande sea \(n\).

Se observa, que las \(n_i\) son las frecuencias observadas de los intervalos \(I_i\) en la muestra, mientras que \(np_{i0}\) son las frecuencias esperadas de los mismos bajo la hipótesis nula, por lo que el estadístico (1.43) es comúnmente referido como

\[ \begin{equation} \chi^{2}=\sum_{i=1}^{k} \frac{\left( O_{i} -E_{i} \right)^2}{E_{i}} \overset{cs}{\rightarrow}\chi_{k-1}^{2}. \tag{1.44} \end{equation} \]

En base al estadístico de prueba \(\chi^{2}\) dado en la ecuación (1.44), la hipótesis alternativa, en este caso, puede ser rechazada si

\[ \begin{equation} \chi^{2}\geq \chi_{k-1;1-\alpha}^{2}. \tag{1.45} \end{equation} \]

Usando el \(p_{valor}\) la hipótesis alternativa es rechazada si la siguiente proposición es verdadera

\[ \begin{equation} p_{valor \ aprox}=1-F_{\chi_{k-1}^{2}}\left(\chi^{2} \right) \leq \alpha. \tag{1.46} \end{equation} \]

Donde \(F_{\chi_{k-1}^{2}}\left(\chi^{2} \right)\) es la función de distribución acumulada de una variable aleatoria ji-cuadrado con \(k-1\) grados de libertad, la cual se encuentra ampliamente tabulada. La distribución base de R proporciona las siguientes funciones, las cuales permiten evaluar estas y otras funciones relacionadas con esta distribución de probabilidad:

dchisq(x, df, ncp = 0, log = FALSE)
pchisq(q, df, ncp = 0, lower.tail = TRUE, log.p = FALSE)
qchisq(p, df, ncp = 0, lower.tail = TRUE, log.p = FALSE)
rchisq(n, df, ncp = 0)

La función dchisq calcula la función de densidad ji-cuadrado, pchisq evalúa la función de distribución ji-cuadrado, qchisq determina los cuantiles de la distribución ji-cuadrado y rchisq genera números aleatorios de esta distribución. Los parámetros de estas funciones se describen en la tabla 1.10:

knitr::kable(
  data.frame(
    stringsAsFactors = FALSE,
    Parámetros = c(
      "x, q", "p", "n", "df",
      "ncp", "log, log.p", "lower.tail"
    ),
    Descripción = c(
      "Vector de cuantiles.",
      "Vector de Probabilidades.",
      "Número de observaciones. Si length(n) > 1, se considera que la longitud es el número requerido.", "Grados de libertad (mayor que 0).",
      "Parámetro de no centralidad (mayor que cero).",
      "Lógico; si es TRUE, las probabilidades p se dan como log(p)",
      "Lógico; si es TRUE (por defecto), las probabilidades son $P[X \\leq x]$, de lo contrario, $P[X > x]$."
    )
  ),
  booktabs = TRUE,
  caption = "Descripción de los parámetros de funciones de R relacionadas con la distribución ji-cuadrado",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2()

Tabla 1.10: Descripción de los parámetros de funciones de R relacionadas con la distribución ji-cuadrado

Parámetros	Descripción
x, q	Vector de cuantiles.
p	Vector de Probabilidades.
n	Número de observaciones. Si length(n) > 1, se considera que la longitud es el número requerido.
df	Grados de libertad (mayor que 0).
ncp	Parámetro de no centralidad (mayor que cero).
log, log.p	Lógico; si es TRUE, las probabilidades p se dan como log(p)
lower.tail	Lógico; si es TRUE (por defecto), las probabilidades son \(P[X \leq x]\), de lo contrario, \(P[X > x]\).

Por otro lado la función chisq.test(x, y = NULL, correct = TRUE, p = rep(1/length(x), length(x)), rescale.p = FALSE, simulate.p.value = FALSE, B = 2000) contenida en la distribución base de R permite ejecutar la prueba de bondad de ajuste ji-cuadrado. Los parámetros de esta función se describen en la tabla 1.11.

knitr::kable(
  data.frame(
    stringsAsFactors = FALSE,
    Parametros = c(
      "x", "y", "correct", "p",
      "rescale.p", "simulate.p.value", "B"
    ),
    Descripcion = c(
      "Un vector o matriz numérica. x e y también pueden ser factores.",
      "Un vector numérico; ignorado si x es una matriz. Si x es un factor, y debe ser un factor de la misma longitud.",
      "Lógica que indica si se debe aplicar o no la corrección de continuidad de Yate al calcular el estadístico de prueba para las tablas de 2 por 2. No se realiza ninguna corrección si simulate.p.value. = TRUE.",
      "Un vector de probabilidades de la misma longitud de x. Se da un error si alguna entrada de p es negativa.",
      "Un parámetro lógico; si es TRUE entonces p es reescalado (si es necesario) para sumar a 1. Si rescale.p es FALSE, y p no suma a 1, se da un error.",
      "Un indicador lógico que indica si se deben calcular los valores p mediante la simulación de Monte Carlo.",
      "Un número entero que especifique el número de réplicas utilizadas en el ensayo de Monte Carlo."
    )
  ),
  booktabs = TRUE,
  caption = "\\label{tab2:par-prueba-jicuadrado}Descripción de los parámetros de la función `chisq.test` de R",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 1.11: Descripción de los parámetros de la función chisq.test de R

Parametros	Descripcion
x	Un vector o matriz numérica. x e y también pueden ser factores.
y	Un vector numérico; ignorado si x es una matriz. Si x es un factor, y debe ser un factor de la misma longitud.
correct	Lógica que indica si se debe aplicar o no la corrección de continuidad de Yate al calcular el estadístico de prueba para las tablas de 2 por 2. No se realiza ninguna corrección si simulate.p.value. = TRUE.
p	Un vector de probabilidades de la misma longitud de x. Se da un error si alguna entrada de p es negativa.
rescale.p	Un parámetro lógico; si es TRUE entonces p es reescalado (si es necesario) para sumar a 1. Si rescale.p es FALSE, y p no suma a 1, se da un error.
simulate.p.value	Un indicador lógico que indica si se deben calcular los valores p mediante la simulación de Monte Carlo.
B	Un número entero que especifique el número de réplicas utilizadas en el ensayo de Monte Carlo.

Una forma más fácil de evaluar la cantidad dada por la ecuación (1.47), es

\[ \begin{equation} \chi ^{2}=\sum_{i=1}^{k}\frac{O_{i}^{2}}{E_i}-n. \tag{1.47} \end{equation} \]

No obstante, es recomendable hacer la determinación a partir del cuadrado de cada diferencia, para detectar cuál es la causante del rechazo, en caso de producirse éste.

Como se dijo al principio de esta sección, se aborda ahora el problema de ajuste a una distribución parcialmente especificada, en la práctica una situación que se presenta con mayor frecuencia. Se empleará la notación \(p_i (\theta_1,\dotsc,\theta_r )\) para las probabilidades de los intervalos \(I_i\) con el fin de resaltar esta diferencia. En uniformidad con el caso anterior, todo se reduce a contrastar, a partir de los datos \(n_i\), con \(i=1,\dotsc,k\), la hipótesis nula \(H_0:\) provienen de una multinomial de parámetros \(n\), \(p_i (\theta_1,\dotsc,\theta_r )\), para \(i=1,\dotsc,k\), con \(k≥r\), frente a la hipótesis alternativa \(H_1:\) provienen de una multinomial de parámetros \(n\), \(p_i (\theta_1,\dotsc,\theta_r )\), para \(i=1,\dotsc,k\).

Si se emplea el contraste de la \(RV\) se debe utilizar

\[ \lambda\left ( n_{1}, \dotsc, n_{k} \right ) =\frac{\underset{\theta_1, \dotsc, \theta_k}{SUP}\frac{n!}{n_1! \dotsm n_k!} p_1^{n_1}\left ( \theta _1, \dotsc,\theta_r \right ) \dotsm p_k^{n_k}\left ( \theta _1, \dotsc,\theta_r \right )}{\underset{p_1, \dotsc, p_k}{SUP}\frac{n!}{n_1! \dotsm n_k!}\left ( p_1 \right )^{n_1} \dotsm \left ( p_k \right )^{n_k}}. \]

Para el cálculo del denominador, se sabe que el supremo se alcanza para \(\widehat{p}_i=\frac{n_i}{n}\), con \(i=1,\dotsc,k\) y para el numerador, hay que resolver en \(θ_1,\dotsc,θ_r\) el sistema

\[ \left.\begin{matrix} \sum_{i=1}^{k}\frac{n_i}{p_i\left ( \theta _1, \dotsc, \theta _r \right )}\frac{\partial }{\partial \theta _j}p_i\left ( \theta _1, \dotsc, \theta _r \right )=0\\ j=1, \dotsc, r \end{matrix}\right\} \]

En el caso de que exista regularidad suficiente, las soluciones de este sistema \(p(\widehat{\theta}_i )\) maximizarán el numerador de la \(RV\), con lo que la región crítica ha de ser

\[ RC= \left\{\prod_{i=1}^{k}\left ( \frac{p_i\left ( \widehat{\theta } \right )}{n_i/n} \right )^{n_i} \leq k \right\} \]

y por los resultados asintóticos de la \(RV\) , bajo \(H_0\), se sigue que

\[ -2\ln\lambda \left ( n_1, \dotsc, n_k \right )\xrightarrow{cs }\chi _{k-r-1}^{2} \]

ya que la hipótesis nula tiene de dimensión \(r\), mientras que \(H_0∪H_1\) tiene de dimensión \(k-1\).

Por lo tanto, la región crítica viene dada por

\[ RC=\left\{\lambda \left ( n_1, \dotsc, n_k \right )=\prod_{i=1}^{k}\left ( \frac{p_i\left ( \widehat{\theta } \right )}{n_i/n} \right )^{n_i}\leq k \right\} \]

Alternativamente, se puede aplicar el método de la distancia de K. Pearson para lo que se introduce la cantidad

\[ D^{2}\left ( p\left ( \widehat{\theta } \right ),\widehat{p} \right )=\sum_{i=1}^{k}\lambda _i\left ( \frac{n_i}{n} -p_i\left ( \widehat{\theta } \right ) \right )^{2} \]

y si se pone \(\lambda_i=\frac{n}{p_i\left ( \widehat{\theta } \right )}\) queda

\[ \begin{equation} D^{2}\left ( p\left ( \widehat{\theta } \right ),\widehat{p} \right )=\sum_{i=1}^{k}\frac{\left ( n_i-np_i\left ( \widehat{\theta } \right ) \right )^{2}}{np_i\left ( \widehat{\theta } \right )}\xrightarrow{cs }\chi _{k-r-1}^{2}. \tag{1.47} \end{equation} \]

Dado que los \(n_i\) son las frecuencias observadas de los intervalos \(I_i\) en la muestra, mientras que \(np_i \left( \widehat{\theta } \right)\) son las frecuencias esperadas estimadas de los mismos bajo la hipótesis nula, es costumbre escribir la ecuación (1.47) como

\[ \begin{equation} \chi^{2} =\sum_{i=1}^{k}\frac{\left ( O_i-\widehat{E}_i \right )^{2}}{\widehat{E}_i}\xrightarrow{cs }\chi _{k-r-1}^{2}. \tag{1.48} \end{equation} \]

Para estimar las frecuencias esperadas de cada clase, es necesario estimar las probabilidades de ocurrencia de cada clase. Para ello es necesario estimar los parámetros de la distribución a la cual se quiere ajustar los datos, a partir de una muestra. La función fitdistr(x, densfun, start, ...) del paquete MASS permite estimar los parámetros de algunas distribuciones univariadas. Los parámetros de la función se describen en la tabla 1.12.

knitr::kable(
  data.frame(
    stringsAsFactors = FALSE,
    Parámetros = c("x", "densfun", "start", "..."),
    Descripción = c(
      "Un vector numérico de longitud mayor o igual a uno que contiene los valores de la muestra.",
      "Se reconocen las distribuciones \"beta\", \"cauchy\", \"chi-squared\", \"exponential\", \"gamma\", \"geometric\", \"lognormal\", \"logistic\", \"negative binomial\", \"normal\", \"Poisson\", \"t\" and \"weibull\", ignorándose otros casos.",
      "Una lista con nombre con los parámetros a optimizar con valores iniciales. Esto puede omitirse para algunas de las distribuciones mencionadas y debe omitirse para otras (véase Detalles).",
      "Parámetros adicionales, tanto para densfun como para optimizar. En particular, se puede utilizar para especificar límites a través de la parte inferior, superior o ambas. Si se incluyen argumentos de densfun (o la función de densidad especificada por la cadena de caracteres correspondientes) se mantendrán fijos."
    )
  ),
  booktabs = TRUE,
  caption = "Descripción de los parámetros de la función fitdistr de R",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 1.12: Descripción de los parámetros de la función fitdistr de R

Parámetros	Descripción
x	Un vector numérico de longitud mayor o igual a uno que contiene los valores de la muestra.
densfun	Se reconocen las distribuciones “beta”, “cauchy”, “chi-squared”, “exponential”, “gamma”, “geometric”, “lognormal”, “logistic”, “negative binomial”, “normal”, “Poisson”, “t” and “weibull”, ignorándose otros casos.
start	Una lista con nombre con los parámetros a optimizar con valores iniciales. Esto puede omitirse para algunas de las distribuciones mencionadas y debe omitirse para otras (véase Detalles).
…	Parámetros adicionales, tanto para densfun como para optimizar. En particular, se puede utilizar para especificar límites a través de la parte inferior, superior o ambas. Si se incluyen argumentos de densfun (o la función de densidad especificada por la cadena de caracteres correspondientes) se mantendrán fijos.

1.4.2 Implementación en R de la Prueba de Bondad de Ajuste Ji-cuadrado

Ejemplo 1.4 (Prueba de Bondad de Ajuste Ji-cuadrado) En 1865 Mendel descubrió un código genético básico mediante la cría de guisantes verdes y amarillos en un experimento. Debido a que el gen amarillo del guisante es dominante, los híbridos de la primera generación fueron amarillos, pero los de la segunda generación eran aproximadamente 75% amarillos y 25% verdes. El color verde reaparece en la segunda generación porque hay un 25% de probabilidad de que dos guisantes, ambos con un gen amarillo y otros con un gen verde, contribuyan con el verde a la siguiente semilla híbrida. En otro experimento de guisantes se consideró tanto el color como las características de la textura. Los resultados de los experimentos se muestran en la tabla 1.13.

datos3 <- data.frame(
  stringsAsFactors = FALSE,
  Tipo.de.Guisante = c(
    "Amarillo Liso", "Amarillo Arrugado",
    "Verde Liso", "Verde Arrugado"
  ),
  Frecuencia.Observada = c(315L, 101L, 108L, 32L),
  Frecuencia.Esperada = c(313L, 104L, 104L, 35L)
)
knitr::kable(
  datos3,
  booktabs = TRUE,
  col.names = c("Tipo de Guisante", "Frecuencia Observada", "Frecuencia Esperada"),
  caption = "\\label{tab2:datos-mendel}Distribución de los guisantes",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2()

Tabla 1.13: Distribución de los guisantes

Tipo de Guisante	Frecuencia Observada	Frecuencia Esperada
Amarillo Liso	315	313
Amarillo Arrugado	101	104
Verde Liso	108	104
Verde Arrugado	32	35

Pruebe la hipótesis de que la ocurrencia de las dos características está determinado por el modelo propuesto (número esperado), con un tamaño \(\alpha=0,01\). Los datos fueron tomados de Kvan y Vidakovic (2007, pág. 157).

En este ejemplo la hipótesis nula que se quiere contrastar es que las frecuencia observados en la muestra se ajustan a las frecuencias esperadas propuestas por Mendel. Ahora, como la distribución de frecuencias propuesta en la hipótesis nula está completamente definida, el estadístico de prueba viene dado por la ecuación (1.44). Luego, evaluando las frecuencias observadas y esperadas en esta ecuación se obtiene el siguiente valor:

\[ \chi^{2}=\sum_{i=1}^{4} \frac{\left( O_{i} -E_{i} \right)^2}{E_{i}} = \frac{\left( 315 -313 \right)^2}{313}+\frac{\left( 101 -104 \right)^2}{104}+\frac{\left( 108 - 104 \right)^2}{104}+\frac{\left( 32 - 35 \right)^2}{35}= 0.510307. \]

Según la ecuación (1.45), el valor del estadístico tabulado es \(\chi_{k-1=3;1-\alpha=0,99}^{2} = 11,3449\) . Mientras que de la ecuación (1.46), el \(p_{valor \, aprox}=1-F_{\chi_{k-1=3}^{2}}\left( 0,5103\right)=1-0,0834 = 0,9166\). Dado que el estadístico de prueba \(\left(0,5103\right)\) es mucho menor que el estadístico tabulado \(\left(13,2767 \right)\) no se puede rechazar la hipótesis de que las frecuencias observadas se ajustan al modelo propuesto por Mendel. A la misma conclusión se llega si se usa el \(p_{valor}\) como criterio de decisión. Los resultados antes descritos se pueden obtener con el siguiente script.

alfa <- 0.01
Ji_cuadrado <- sum((datos3$Frecuencia.Observada - datos3$Frecuencia.Esperada)^2 /
  datos3$Frecuencia.Esperada)
est_tab <- qchisq(p = 1 - alfa, df = nrow(datos3) - 1)
p_valor <- 1 - pchisq(q = Ji_cuadrado, df = nrow(datos3) - 1)
paste("El estadísitico Ji-cuadrado =", round(Ji_cuadrado, digits = 4))
paste("El estadístico tabulado =", round(est_tab, digits = 4))
paste("El p_valor =", round(p_valor, digits = 4))

#> [1] "El estadísitico Ji-cuadrado = 0.5103"
#> [1] "El estadístico tabulado = 11.3449"
#> [1] "El p_valor = 0.9166"

El estadístico de prueba \(\chi^2\) y el \(p_{valor \, aprox.}\) se pueden obtener de manera directa con la función chisq.test descrita en la tabla 1.11, como se muestra en el siguiente código.

chisq.test(
  x = datos3$Frecuencia.Observada,
  p = datos3$Frecuencia.Esperada /
    sum(datos3$Frecuencia.Esperada)
)

#> 
#>  Chi-squared test for given probabilities
#> 
#> data:  datos3$Frecuencia.Observada
#> X-squared = 0.51031, df = 3, p-value = 0.9166

Ejemplo 1.5 (Prueba de Bondad de Ajuste Ji-cuadrado) Se ha estimado que el número de accidentes diarios en cada regimiento del ejército sigue una distribución de Poisson de parámetro 2. Un determinado regimiento ha recogido, durante 200 días, los siguientes datos:

\[ \begin{array}{lccccccc} \hline N^{\circ} \, de \, accidentes & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 \\ N^{\circ} \, de \, días & 22 & 53 & 58 & 39 & 20 & 5 & 2 & 1 \\ \hline \end{array} \]
Los datos fueron tomados de Vélez y García (2012, pág. 435).

En la tabla 1.14 se aprecian las frecuencia observadas y esperadas. Donde las frecuencias esperadas viene dada por \(E_i=np_{0i}\). Donde \(n\) es el tamaño de la muestra (200) y \(p_{0i}\) son las probabilidades esperadas, las cuales se obtiene sustituyendo en la función de probabilidad de Poisson

\[ f\left(x, \lambda \right)=e^{\lambda}\frac{\lambda^x}{x!}, \]

\(x=0, 1, 2, \dotsc7\) y \(\lambda=2\). La distribución base de R tiene la función dpois que permite evaluar esta función.

datos4 <- tibble(
  accidentes = c(0L, 1L, 2L, 3L, 4L, 5L, 6L, 7L),
  dias = c(22L, 53L, 58L, 39L, 20L, 5L, 2L, 1L)
) %>%
  mutate(
    p0i = dpois(x = accidentes, lambda = 2),
    Ei = sum(dias) * p0i
  )
knitr::kable(
  datos4,
  digits = 4,
  booktabs = TRUE,
  row.names = TRUE,
  col.names = c("$N^{\\circ} \\, de \\, accidentes$", "$N^{\\circ} \\, de \\, días \\left( O_i \\right)$", "$p_{0i}$", "$E_i$"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  caption = "\\label{tab2:ejemplo4}Frecuencias obsevadas y esperadas",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  row_spec(which(datos4$Ei < 5), color = "red", strikeout = T) %>%
  kable_classic_2() 

Tabla 1.14: Frecuencias obsevadas y esperadas

	\(N^{\circ} \, de \, accidentes\)	\(N^{\circ} \, de \, días \left( O_i \right)\)	\(p_{0i}\)	\(E_i\)
1	0	22	0,1353	27,0671
2	1	53	0,2707	54,1341
3	2	58	0,2707	54,1341
4	3	39	0,1804	36,0894
5	4	20	0,0902	18,0447
6	5	5	0,0361	7,2179
7	6	2	0,0120	2,4060
8	7	1	0,0034	0,6874

Como se nota en la tabla anterior, las clases 7 y 8 tienen frecuencias esperadas menores que 5, y peor aún, la de la clase 8 es menor que 1. Por tanto, estas se funden con la clase 6 para dar lugar a la clase 5 o más accidentes en un día, para que sea procedente utilizar la prueba de bondad de ajuste Ji-Cuadrado.

datos4 <- bind_rows(datos4[1:5, ], colSums(datos4[6:8, ])) %>%
  mutate(accidentes = as.character(accidentes))
datos4[6, 1] <- "$\\geq5$"
knitr::kable(
  datos4,
  digits = 4,
  booktabs = TRUE,
  row.names = TRUE,
  col.names = c("$N^{\\circ} \\, de \\, accidentes$", "$N^{\\circ} \\, de \\, días \\left( O_i \\right)$", "$p_{0i}$", "$E_i$"),
  align = c("r"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  caption = NULL,
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2() 

	\(N^{\circ} \, de \, accidentes\)	\(N^{\circ} \, de \, días \left( O_i \right)\)	\(p_{0i}\)	\(E_i\)
1	0	22	0,1353	27,0671
2	1	53	0,2707	54,1341
3	2	58	0,2707	54,1341
4	3	39	0,1804	36,0894
5	4	20	0,0902	18,0447
6	\(\geq5\)	8	0,0516	10,3113

Dado que con este ajuste se ha logrado que todas las frecuencias esperadas sean mayor que 5 y como la distribución que se pretende ajustar a los datos está completamente especificada. Entonces el estadístico calculado, según la ecuación (1.44), viene dado por:

\[ \chi^{2}=\sum_{i=1}^{4} \frac{\left( O_{i} -E_{i} \right)^2}{E_{i}} = \frac{\left( 22 - 27,06 \right)^2}{27,06}+\frac{\left(53 -54,14 \right)^2}{54,14}+\frac{\left(58 - 54,14 \right)^2}{54,14}+\frac{\left(39 - 36,08 \right)^2}{36,08}+\frac{\left(20 - 18,04 \right)^2}{18,04}+\frac{\left(8 - 10,30 \right)^2}{10,30}= 2.2130844. \]

De manera análoga al ejemplo 3, de la ecuación (1.46), el \(p_{valor \, aprox}\) es:

\[ p_{valor \, aprox}=1-F_{\chi_{k-1=5}^{2}}\left( 2,2131 \right)=1- 0,1811 = 0,8189. \]

De los resultados anteriores, se concluye que no se puede rechazar la hipótesis de que los accidentes diarios en cada regimiento del ejercito sigue una distribución de Poisson con \(\lambda = 2\), a un nivel de significancia del 5%.

chisq.test(
  x = datos4$dias, p = datos4$Ei,
  rescale.p = TRUE
)

#> 
#>  Chi-squared test for given probabilities
#> 
#> data:  datos4$dias
#> X-squared = 2.2104, df = 5, p-value = 0.8193

Ejemplo 1.6 (Prueba de Bondad de Ajuste Ji-cuadrado) Una máquina, en correcto estado de funcionamiento, produce piezas cuya longitud se distribuye normalmente con media 10,5 y desviación estándar 0,15. En determinado momento se observa la siguiente muestra, de tamaño 40, de la longitud de las piezas producidas:

\[ \begin{array}{ccccccc} \hline 10,39 & 10,66 & 10,12 & 10,32 & 10,25 & 10,91 & 10,52 & 10,83\\ 10,72 & 10,28 & 10,35 & 10,46 & 10,54 & 10,72 & 10,23 & 10,18\\ 10,62 & 10,49 & 10,32 & 10,61 & 10,64 & 10,23 & 10,29 & 10,78\\ 10,81 & 10,39 & 10,34 & 10,62 & 10,75 & 10,34 & 10,41 & 10,81\\ 10,64 & 10,53 & 10,31 & 10,46 & 10,47 & 10,43 & 10,57 & 10,74\\ \hline \end{array} \]

y se desea saber si la muestra avala que la máquina está funcionando correctamente. Los datos fueron tomados de Vélez y García (2012, pág. 437).

Para aplicar el contraste de bondad de ajuste ji-cuadrado, hay que construir la tabla de distribución de frecuencia en intervalos de clases; como \(n=40\), se pueden hacer 8 intervalos, con números esperado de observaciones igual a 5, buscando los cuantiles de órdenes \(0,125, 0,25, 0,375, \dotsc, 0,875\). El siguiente código permite ejecutar el procedimiento descrito previamente.

lp <- c(
  10.39, 10.66, 10.12, 10.32, 10.25, 10.91, 10.52, 10.83,
  10.72, 10.28, 10.35, 10.46, 10.54, 10.72, 10.23, 10.18,
  10.62, 10.49, 10.32, 10.61, 10.64, 10.23, 10.29, 10.78,
  10.81, 10.39, 10.34, 10.62, 10.75, 10.34, 10.41, 10.81,
  10.64, 10.53, 10.31, 10.46, 10.47, 10.43, 10.57, 10.74
)
datos5 <- as.data.frame(table(cut(
  x = lp, breaks =
    qnorm(
      p = seq(
        from = 0, to = 1,
        by = 0.125
      ),
      mean = 10.5, sd = 0.15
    )
)))
datos5 <- datos5 %>%
  transmute(
    Clase = str_replace_all(
      string = as.character(Var1),
      c(
        "\\(" = "$\\\\left(",
        "," = ";\\\\:",
        "\\." = ",",
        "Inf" = "\\\\infty",
        "\\]" = "\\\\right]$"
      )
    ),
    Oi = Freq,
    p0i = rep(x = 1 / 8, times = 8),
    Ei = p0i * sum(Oi)
  )
knitr::kable(
  datos5,
  booktabs = TRUE,
  row.names = TRUE,
  col.names = c("Clase", "$O_i$", "$p_{0i}$", "$E_i$"),
  align = c("r"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  caption = NULL,
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2()

	Clase	\(O_i\)	\(p_{0i}\)	\(E_i\)
1	\(\left(-\infty;\:10,33\right]\)	10	0,125	5
2	\(\left(10,33;\:10,4\right]\)	5	0,125	5
3	\(\left(10,4;\:10,45\right]\)	2	0,125	5
4	\(\left(10,45;\:10,5\right]\)	4	0,125	5
5	\(\left(10,5;\:10,55\right]\)	3	0,125	5
6	\(\left(10,55;\:10,6\right]\)	1	0,125	5
7	\(\left(10,6;\:10,67\right]\)	6	0,125	5
8	\(\left(10,67;\: \infty\right]\)	9	0,125	5

El siguiente código muestra la ejecución de la prueba de bondad de ajuste Ji-cuadrado para contrastar la hipótesis de que la longitud de la pieza se distribuye normalmente con media 10,5 y desviación estándar 0,15. Note que no se ha especificado el valor del parámetro p de la función chisq.test, por lo tanto este asume el valor por defecto descrito en la tabla 1.11.

Como se observa en la salida, el \(p_{valor\, aprox} = 0.0445075\) lo que indica que se rechaza la hipótesis de que la longitud de la pieza se ajusta a una distribución normal con media 10,5 y desviación estándar 0,15; con un nivel de confianza del 95%.

chisq.test(x = datos5$Oi)

#> 
#>  Chi-squared test for given probabilities
#> 
#> data:  datos5$Oi
#> X-squared = 14.4, df = 7, p-value = 0.04451

La siguiente gráfica corrobora la conclusión de que los datos no se ajustan a una distribución \(N\left(\mu = 10,5;\: \sigma = 0,15 \right)\). Como se ve en la figura 1.1 la densidad \(N\left(\mu = 10,5; \sigma = 0,15 \right)\) , representada por la línea roja, no se ajusta a la densidad empírica de los datos muestrales, representada por la línea azul.

ggplot(as.data.frame(lp), aes(x = lp)) +
  geom_histogram(aes(y = ..density..),
    bins = 8, color = "white", fill = "steelblue"
  ) +
  stat_function(
    fun = dnorm, color = "red",
    args = list(mean = 10.5, sd = 0.15)
  ) +
  geom_line(stat = "density", color = "blue") +
  expand_limits(y = 0) +
  xlab("Longitud de la pieza") +
  ylab("Densidad") +
  theme_minimal()

Figura 1.1: Densidad de la longitud de la pieza

<br>

Ejemplo 1.7 (Prueba de Bondad de Ajuste Ji-cuadrado) La muestra del ejemplo 1.6 ha puesto de relieve un funcionamiento incorrecto de la máquina, debido a que la longitud de las piezas producidas por esta no se ajustan a la distribución \(N\left(\mu = 10,5; \: \sigma = 0,15 \right)\). Cabe, sin embargo, la posibilidad de que el desajuste haya afectado a la media y la desviación típica, pero que la distribución de la longitud de las piezas producidas siga siendo de tipo normal. Para verificar tal posibilidad hay que estimar (por máxima verosimilitud) la media y la desviación típica poblacional.

El estimación por máxima verosimilitud de la media de la distribución normal es:

\[ \hat{\mu} = \bar{x} = \frac{\sum_{i=1}^{n}xi}{n} = \frac{420,08}{40} = 10,502. \]

Mientras que la estimación de la desviación estándar por este método viene dada por:

\[ \hat{\sigma }=\sqrt{\frac{\sum_{i=1}^{n}\left ( x_i - \bar{x} \right )^{2}}{n}}=\sqrt{\frac{1,64144}{40}}= 0,2025734. \]

La función fitdistr del paquete MASS cuyos parámetro se describen en la tabla 1.12 calcula los estimadores de máxima verosimilitud de algunas distribuciones de probabilidad univariadas, entre ellas la distribución normal. Con el siguiente código se determinan tales estimaciones, a partir de la información de la muestra.

fitdistr(x = lp, densfun = "normal")

#>       mean           sd     
#>   10.50200000    0.20257344 
#>  ( 0.03202967) ( 0.02264840)

El siguiente código permite construir la tabla de distribución de frecuencia observadas y esperadas estimadas. Estas últimas se calculan a partir de las probabilidades estimadas de cada clase (\(\hat{E}_i=n\hat{p}_{0i}\)). La construcción de la tabla se obtiene con el siguiente código.

lim_int <- c(-Inf, 10.2, 10.3, 10.4, 10.5, 10.6, 10.7, 
             10.8, 10.9, Inf)
datos6 <- as.data.frame(table(cut(
  x = lp, breaks = lim_int
)))
datos6 <- datos6 %>%
  transmute(
    Clase = str_replace_all(
      string = as.character(Var1),
      c(
        "\\(" = "$\\\\left(",
        "," = ";\\\\:",
        "\\." = ",",
        "Inf" = "\\\\infty",
        "\\]" = "\\\\right]$"
      )
    ),
    Oi = Freq,
    p0i = diff(pnorm(q = lim_int, mean = 10.502, sd = 0.2026)),
    Ei = p0i * sum(Oi)
  )
knitr::kable(
  datos6,
  digits = 4,
  booktabs = TRUE,
  row.names = TRUE,
  col.names = c("Clase", "$O_i$", "$\\hat{p}_{0i}$", "$\\hat{E}_i$"),
  align = c("r"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  caption = NULL,
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  row_spec(which(datos6$Ei < 5), color = "red", strikeout = T) %>%
  kable_classic_2()

	Clase	\(O_i\)	\(\hat{p}_{0i}\)	\(\hat{E}_i\)
1	\(\left(-\infty;\:10,2\right]\)	2	0,0680	2,7212
2	\(\left(10,2;\:10,3\right]\)	5	0,0913	3,6537
3	\(\left(10,3;\:10,4\right]\)	8	0,1479	5,9180
4	\(\left(10,4;\:10,5\right]\)	6	0,1887	7,5496
5	\(\left(10,5;\:10,6\right]\)	4	0,1896	7,5857
6	\(\left(10,6;\:10,7\right]\)	6	0,1501	6,0033
7	\(\left(10,7;\:10,8\right]\)	5	0,0935	3,7420
8	\(\left(10,8;\:10,9\right]\)	3	0,0459	1,8369
9	\(\left(10,9;\: \infty\right]\)	1	0,0247	0,9895

Como se puede ver en la tabla anterior, las dos primeras clases tienen frecuencias esperadas estimadas menores que 5, por la que estas se funden en una sola para dar lugar a la clase longitud de la pieza menor o igual a 10,3. Del mismo modo, las tres últimas clases tienen frecuencias esperadas estimadas menores que 5 por lo se colapsan para dar lugar a la clase longitud de la pieza mayor que 10,7. La nueva tabla se muestra con el siguiente código.

lim_int <- c(-Inf, 10.3, 10.4, 10.5, 10.6, 10.7, Inf)
datos6 <- as.data.frame(table(cut(
  x = lp, breaks = lim_int
)))
datos6 <- datos6 %>%
  transmute(
    Clase = str_replace_all(
      string = as.character(Var1),
      c(
        "\\(" = "$\\\\left(",
        "," = ";\\\\:",
        "\\." = ",",
        "Inf" = "\\\\infty",
        "\\]" = "\\\\right]$"
      )
    ),
    Oi = Freq,
    p0i = diff(pnorm(q = lim_int, mean = 10.502, sd = 0.2026)),
    Ei = p0i * sum(Oi)
  )
knitr::kable(
  datos6,
  digits = 4,
  booktabs = TRUE,
  row.names = TRUE,
  col.names = c("Clase", "$O_i$", "$\\hat{p}_{0i}$", "$\\hat{E}_i$"),
  align = c("r"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  caption = NULL,
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2()

	Clase	\(O_i\)	\(\hat{p}_{0i}\)	\(\hat{E}_i\)
1	\(\left(-\infty;\:10,3\right]\)	7	0,1594	6,3749
2	\(\left(10,3;\:10,4\right]\)	8	0,1479	5,9180
3	\(\left(10,4;\:10,5\right]\)	6	0,1887	7,5496
4	\(\left(10,5;\:10,6\right]\)	4	0,1896	7,5857
5	\(\left(10,6;\:10,7\right]\)	6	0,1501	6,0033
6	\(\left(10,7;\: \infty\right]\)	9	0,1642	6,5685

Fíjese que la prueba de hipótesis que se plantea en este ejemplo es compuesta, dado que los parámetros de la distribución bajo la hipótesis nula son desconocidos (la media \(\mu\) y la desviación estándar \(\sigma\)). Por lo tanto, el estadístico de la prueba se distribuye asintóticamente Ji-cuadrado con \(k-r-1=6-2-1=3\) grados de libertad. El valor del estadístico \(\chi^{2}\) y el \(p_{valor \, aprox}\) se generan con el siguiente código.

X2 <- sum((datos6$Oi - datos6$Ei)^2 / datos6$Ei)
p.valor <- pchisq(q = X2, df = nrow(datos6) - 2 - 1, lower.tail = F)
paste("El estadístic Ji-Cuadrado =", X2)
paste("el p.valor =", p.valor)

#> [1] "El estadístic Ji-Cuadrado = 3.70690303455726"
#> [1] "el p.valor = 0.294902155121217"

Dado que el \(p_{valor \, aprox} = 0,2949022\) es mayor que 0,05 no se puede rechazar la hipótesis de que los datos proviene de una población normal.

Si se propone ahora realizar la prueba, de tal manera que la tabla de frecuencia no amerite ser corregida, como se realizó en el ejemplo 5. La tabla quedaría de la siguiente manera, ejecutando el siguiente script.

datos6 <- as.data.frame(table(cut(
  x = lp, breaks =
    qnorm(
      p = seq(
        from = 0, to = 1,
        by = 0.125
      ),
      mean = 10.502, sd = 0.2026
    )
)))
datos6 <- datos6 %>%
  transmute(
    Clase = str_replace_all(
      string = as.character(Var1),
      c(
        "\\(" = "$\\\\left(",
        "," = ";\\\\:",
        "\\." = ",",
        "Inf" = "\\\\infty",
        "\\]" = "\\\\right]$"
      )
    ),
    Oi = Freq,
    p0i = rep(x = 1 / 8, times = 8),
    Ei = p0i * sum(Oi)
  )
knitr::kable(
  datos6,
  booktabs = TRUE,
  row.names = TRUE,
  col.names = c("Clase", "$O_i$", "$p_{0i}$", "$E_i$"),
  align = c("r"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  caption = NULL,
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2()

	Clase	\(O_i\)	\(p_{0i}\)	\(E_i\)
1	\(\left(-\infty;\:10,27\right]\)	5	0,125	5
2	\(\left(10,27;\:10,37\right]\)	8	0,125	5
3	\(\left(10,37;\:10,44\right]\)	4	0,125	5
4	\(\left(10,44;\:10,5\right]\)	4	0,125	5
5	\(\left(10,5;\:10,57\right]\)	3	0,125	5
6	\(\left(10,57;\:10,64\right]\)	4	0,125	5
7	\(\left(10,64;\:10,74\right]\)	5	0,125	5
8	\(\left(10,74;\: \infty\right]\)	7	0,125	5

Como se puede ver en la tabla anterior, todas las frecuencias esperadas son mayores o iguales que 5, por lo que se procede a aplicar la prueba, por medio de las siguientes instrucciones.

X2 <- sum((datos6$Oi - datos6$Ei)^2 / datos6$Ei)
p.valor <- pchisq(q = X2, df = nrow(datos6) - 2 - 1, lower.tail = F)
paste("El estadístic Ji-Cuadrado =", X2)
paste("el p.valor =", p.valor)

#> [1] "El estadístic Ji-Cuadrado = 4"
#> [1] "el p.valor = 0.54941595135278"

Como se nota en la salida anterior el \(p_{valor \, aprox} = 0,549416\) lo que conduce, de igual manera, a aceptar la hipótesis de que los datos se distribuyen normalmente.

1.4.3 Problemas Propuestos

1. Juan Pérez, vendedor de la empresa La Margariteña, S.A, tiene siete empresas que visita semanalmente. Se piensa que las ventas del señor Juan pueden describirse mediante la distribución binomial, con probabilidad de venta de 0,45. Examinando la distribución de frecuencia observada del número de ventas por semana del señor Juan, determine si la distribución corresponde en efecto a la distribución sugerida. Use el nivel de significancia de 5%.

\[ \begin{array}{lcccccccc} \hline Número\,de\,ventas\,por\,semana & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7\\ Frecuencia\,del\,número\,de\,ventas & 25 & & 32 & 61 & 47 & 39 & 21 & 18 & 12\\ \hline \end{array} \]

2. Los vigilantes que monitorean una entrada de seguridad a una empresa están teniendo problemas con las largas filas que se forman durante las horas de tráfico pico. Antes de continuar su análisis del problema, necesitan saber la distribución aproximada de la llegada de carros. Recabaron los siguientes datos cada minuto entre las 7:10 y las 8:00.

\[ \begin{array}{lccccccccccc} \hline Número\, de\, llegadas & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11\\ Frecuencia & 5 & 3 & 2 & 6 & 6 & 2 & 6 & 10 & 4 & 4 & 2\\ \hline \end{array} \]

Pruebe si una distribución de Poisson con una media de seis describe adecuadamente estos datos, use el nivel de significancia de 5%.

3. Los impactos de 60 bombas volantes lanzadas por los alemanes, en la segunda guerra mundial, sobre la superficie de Londres, considerada cuadrada, fueron clasificados en 9 zonas obtenidas dividiendo cada lado en tres partes iguales, con los siguientes resultados:

\[ \begin{array}{ccc} \hline 8 & 7 & 3\\ 5 & 9 & 11\\ 6 & 4 & 7\\ \hline \end{array} \]

Los responsables de la defensa querían averiguar si las bombas tenían algún objetivo concreto o se distribuían al azar sobre la superficie de la ciudad. Los datos fueron tomados de Vélez y García (2012, pág. 438).

1.5 Prueba de Bondad de Ajuste Kolgomorov-Smirnov

La prueba de bondad de ajuste Kolgomorov-Smirnov se usa para contrastar la hipótesis de que una variable continua tiene cierta distribución de probabilidad. Como se mostró en la sección 1.4, la prueba Ji-cuadrado discretiza las observaciones muestrales en conjuntos de una cierta partición, para comparar después el histograma de frecuencias observadas con el histograma de frecuencias que debería obtenerse de acuerdo con la hipótesis a contrastar. De esta manera, en el caso de distribuciones poblacionales continuas, dicho test no hace el mejor uso posible de la información contenida en la muestra, puesto que ignora el valor exacto de las observaciones. Con un tamaño muestral relativamente grande, que permita utilizar una partición fina del conjunto de observaciones, el efecto puede no ser acusado. Pero con un número pequeño de observaciones y, por tanto, con una partición grosera de ellas, la pérdida de información puede ser importante; además, en tal caso, la distribución exacta del estadístico de contraste puede estar alejada de la distribución \(\chi^{2}\) límite.

Así pues, en el contraste de bondad de ajuste de una distribución teórica unidimensional de tipo continuo, es a menudo preferible el uso de la prueba Kolgomorov-Smirnov.

1.5.1 Racionalización de la Prueba de Bondad de Ajuste Kolgomorov-Smirnov

Sea \(x_1,\dotsc,x_n\) una realización particular de una m.a.s. de tamaño \(n\) de una población \(X\) con función de distribución \(F\left(x \right)\) continua, la función de distribución empírica o muestral se define mediante

\[ \begin{equation} F_{n}^{*}\left ( x \right )=\begin{cases} 0, & \text{ si } x< x_{\left ( i \right )}\\ \frac{i}{n}, & \text{ si } x_{\left ( i \right )}\leq x< x_{\left ( i +1\right )}, \text{ para }\: i=1, \dotsc, n-1 \\ 1, & \text{ si } x\geq x_{\left ( n \right )} \end{cases}. \tag{1.49} \end{equation} \]

Donde \(x_{\left ( 1 \right )},\dotsc,x_{\left ( n \right )}\) es la muestra ordenada de menor a mayor, con \(i= 0,1,\dotsc,n\); donde aquí se añade por convenio \(x_{\left ( 0 \right )}=-\infty\) y \(x_{\left ( n+1 \right )}=-\infty\). Alternativamente la función \(F_{n}^{*}\left ( x \right )\) puede escribirse como:

\[ \begin{equation} F_{n}^{*}\left ( x \right )=\frac{1}{n}\sum_{i=1}^{n}I_{(-\infty , x]}\left ( x_{i} \right ) \tag{1.50} \end{equation} \]

La función ecdf determina la función de distribución acumulada empírica de una muestra dada por la ecuación (1.49), mientras que la función Ecdf del paquete Hmisc la calcula y la grafica.

Definición 1.1 (Estadísticos de Kolmogorov-Smirnov) Se llaman estadísticos unilaterales de Kolmogorov-Smirnov a

\[ \begin{equation} \begin{split} D_{n}^{+} & = \underset{x}{Sup}\left ( F_{n}^{*}\left ( x \right )-F\left ( x \right ) \right )\\ D_{n}^{-} & = \underset{x}{Sup}\left ( F\left ( x \right ) - F_{n}^{*} (x)\right ) \end{split} \tag{1.51} \end{equation} \]

y estadístico bilateral de Kolmogorov-Smirnov a

\[ \begin{equation} D_{n}=\underset{x}{Sup}\left| F_{n}^{*}\left ( x \right ) -F\left ( x \right )\right| \tag{1.52} \end{equation} \]

Para demostrar que los estadísticos de Kolmogorov-Smirnov son de distribución libre se hace uso del siguiente teorema.

Teorema 1.1 (Glivenko-Cantelli) Si se tiene una m.a.s. de tamaño \(n\) de una población \(X\), con función de distribución \(F\left ( x \right )\) , para cualquier número real positivo arbitrario \(\varepsilon\), se tiene que

\[ \begin{equation} \displaystyle \lim_{n \to \infty}P\left\{ \underset{x\in \mathbb{R}}{Sup}\left| F_{n}^{*}\left ( x \right )-F\left ( x \right )\right|\geq \varepsilon \right\}=0. \tag{1.53} \end{equation} \]

El teorema de Glivenko-Cantelli (teorema 1.1), asegura que \(F_{n}^{*}\left ( x \right )\) converge a \(F\left ( x \right )\) uniformemente con probabilidad uno. Por lo tanto la magnitud de las diferencias entre \(F_{n}^{*}\left ( x \right )\) y \(F\left ( x \right )\), proporciona información de la compatibilidad entre la muestra y \(F\left ( x \right )\).

Proposición 1.1 Para una m.a.s. de tamaño \(n\) de una población \(X\) con función de distribución \(F\left ( x \right )\) continua, la función de distribución de los estadísticos \(D_{n}^{+}\), \(D_{n}^{-}\) y \(D_{n}\) no depende de \(F\).

Demostración. Recordando que \(x_{\left ( 0 \right )}=-\infty\) y \(x_{\left ( n+1 \right )}=-\infty\), se puede poner

\[ D_{n}^{+}=\underset{0\leq i\leq n}{m\acute{a}x}\:\:\underset{x_{\left ( i \right )}\leq x\leq x_{\left ( i+1 \right )}}{Sup}\left ( \frac{i}{n} -F\left ( x \right )\right ) \]

Como en el intervalo \(\left [ x_{\left ( n \right )},\: x_{\left ( n+1 \right )}\right )\) la función \(1-F\left ( x \right )\) es positiva o nula y \(\underset{x_{\left ( i \right )}\leq x< x_{\left ( i+1 \right )}}{Sup}\left ( \frac{i}{n} -F\left ( x_{\left ( i \right )} \right )\right )=\frac{i}{n}-F\left ( x_{\left ( i \right )} \right )\) se tiene que

\[ \begin{equation} D_{n}^{+}=\underset{1\leq i\leq n}{m\acute{a}x}\left ( \frac{i}{n} -F\left ( x_{\left ( i \right )} \right )\right ) \tag{1.54} \end{equation} \]

como \(F\left ( x_{\left ( i \right )} \right )=U_{\left ( i \right )}\) es el estadístico ordenado de orden \(i\) de una muestra de la v.a. \(U\left ( a=0, b=1 \right )\); es claro que esta expresión no depende de la distribución \(F\).

Análogamente

\[ \begin{split} D_{n}^{-}&=\underset{0\leq i\leq n}{m\acute{a}x} \:\:\underset{x_{\left ( i \right )}\leq x\leq x_{\left ( i+1 \right )}}{Sup}\left ( F\left ( x \right ) -\frac{i}{n}\right )\\ &=\underset{0\leq i\leq n}{m\acute{a}x}\left ( F\left ( x_{\left ( i+1 \right )} \right )-\frac{i}{n} \right )\\ &=\underset{1\leq i\leq n}{m\acute{a}x}\left ( F\left ( x_{\left ( i \right )} \right )-\frac{i-1}{n} \right )\\ \end{split} \]

sin más que observar que en el intervalo \(\left[x_{\left ( 0 \right )},x_{\left ( 1 \right )}\right)\) la diferencia \(F\left ( x_{\left ( i+1 \right )} \right )-\frac{i}{n}\) es positiva o nula.

En Resumen,

\[ \begin{equation} D_{n}^{-}=\underset{1\leq i\leq n}{m\acute{a}x}\left ( F\left ( x_{\left ( i \right )} \right )-\frac{i-1}{n} \right ). \tag{1.55} \end{equation} \]

Por último al ser

\[ \begin{equation} D_{n}=m\acute{a}x\left\{D_{n}^{+}, D_{n}^{-}\right\} \tag{1.56} \end{equation} \]

se sigue que la distribución de los tres estadísticos es independiente de la función de distribución \(F\).

Se puede entonces enunciar esta proposición diciendo, que los estadísticos de Kolmogorov-Smirnov tienen distribución libre.

Las distribuciones asintóticas de estos estadísticos fueron obtenidas por Kolmogorov (1933, 1941) y por Smirnov (1939, 1948), quien dió una demostración más sencilla y generalizó el problema al caso de dos muestras.

La distribución de probabilidad de los estadísticos de Kolgomorov-Smirnov se encuentra ampliamente tabulada para diferentes valores de \(n\) y \(\alpha\). En particular, las funciones pkolm y pkolmin del paquete kolmim permite evaluar la distribución de Kolgomorov.

La utilización de estos estadísticos, para contrastar si una m.a.s. proviene de una distribución totalmente especificada \(F_0\left(x\right)\), queda recogida en la tabla 1.15.

knitr::kable(
  data.frame(
    stringsAsFactors = FALSE,
    H0 = c("", "$F\\left(x \\right)=F_0\\left(x \\right)$", ""),
    H1 = c(
      "$F\\left(x \\right)\\neq F_0\\left(x \\right)$",
      "$F\\left(x \\right)>F_0\\left(x \\right)$",
      "$F\\left(x \\right)<F_0\\left(x \\right)$"
    ),
    region_critica = c(
      "$D_{n}> d_{n,\\alpha }$",
      "$D_{n}^{+}> d_{n,\\alpha}^{+}$",
      "$D_{n}^{-}> d_{n,\\alpha }^{-}$"
    ),
    p_valor = c(
      "$P \\left( D_{n} > d_{n} \\right)$",
      "$P \\left(D_{n}^{+ }> d_{n}^{+} \\right)$",
      "$P \\left(D_{n}^{-}> d_{n}^{-} \\right)$"
    )
  ),
  booktabs = TRUE,
  row.names = FALSE,
  col.names = c("$H_0$", "$H_1$", "Región Crítica", "$p_{valor}$"),
  align = c("c", "c", "c"),
  caption = "\\label{tab2:regioncritica-kolgomorovsmirnov}Regiones críticas del test de Kolmogorov-Smirnov para una muestra",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2()

Tabla 1.15: Regiones críticas del test de Kolmogorov-Smirnov para una muestra

\(H_0\)	\(H_1\)	Región Crítica	\(p_{valor}\)
	\(F\left(x \right)\neq F_0\left(x \right)\)	\(D_{n}> d_{n,\alpha }\)	\(P \left( D_{n} > d_{n} \right)\)
\(F\left(x \right)=F_0\left(x \right)\)	\(F\left(x \right)>F_0\left(x \right)\)	\(D_{n}^{+}> d_{n,\alpha}^{+}\)	\(P \left(D_{n}^{+ }> d_{n}^{+} \right)\)
	\(F\left(x \right)<F_0\left(x \right)\)	\(D_{n}^{-}> d_{n,\alpha }^{-}\)	\(P \left(D_{n}^{-}> d_{n}^{-} \right)\)

La función ks.test(x, y, …, alternative = c("two.sided", "less", "greater"), exact = NULL) de la distribución base de R realiza la prueba de bondad de ajuste Kolmogorov-Smirnov de una o dos muestras. Los parámetros de esta función se describen en la tabla 1.16.

knitr::kable(
  data.frame(
    stringsAsFactors = FALSE,
    Parámetros = c("x", "y", "…", "alternative"),
    Descripción = c(
      "Un vector numérico de datos.",
      "Un vector numérico de datos o una cadena de caracteres que designa una función de distribución acumulada o una función de distribución acumulada conocida como por ejemplo pnorm. Sólo son válidas las FDA continuas",
      "Parámetros de la distribución especificada en y (como una cadena de caracteres)",
      "Indica la hipótesis alternativa y debe ser: \" \"two.sided\" (por defecto), \" less \" o \" greater \". Puede especificar sólo la letra inicial del valor, pero el nombre del argumento debe estar completo"
    )
  ),
  booktabs = TRUE,
  caption = "Descripción de los parámetros de la función ks.test de R",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2()

Tabla 1.16: Descripción de los parámetros de la función ks.test de R

Parámetros	Descripción
x	Un vector numérico de datos.
y	Un vector numérico de datos o una cadena de caracteres que designa una función de distribución acumulada o una función de distribución acumulada conocida como por ejemplo pnorm. Sólo son válidas las FDA continuas
…	Parámetros de la distribución especificada en y (como una cadena de caracteres)
alternative	Indica la hipótesis alternativa y debe ser: " “two.sided” (por defecto), " less " o " greater ". Puede especificar sólo la letra inicial del valor, pero el nombre del argumento debe estar completo

De igual manera, la función ks.test.imp(x, y, ...) del paquete kolmin permite ejecutar esta prueba.

1.5.2 Implementación en R de la Prueba de Bondad de Ajuste Kolmogorov-Smirnov

Ejemplo 1.8 (Prueba de Bondad de Ajuste Kolmogorov-Smirnov) Consideremos un experimento industrial en el que tenemos una muestra de n=10 tejidos sometido a un lavado. El objetivo del experimento es analizar el desempeño de un nuevo detergente experimental para la ropa. Específicamente, la variable de respuesta bajo estudio es la llamada reflectancia, es decir, la proporción de luz incidente que una determinada superficie (de tela) es capaz de reflejar, lo que puede ser considerado una medida relacionada con la eficacia de limpieza del detergente. Supongamos que deseamos probar, en el nivel de significación \(\alpha=5\%\), si la reflectancia está o no uniformemente distribuida, es decir, si sigue la ley de distribución U(0,1). Los datos fueron tomados de Bonnini et al (2014, pág. 28).

\[ \begin{array}{ccccccc} \hline \mathbf{Pieza \, de\, Tela} & 1 & 2 & 3 & 4 & 5\\ \mathbf{Reflectancia} & 0,608 & 0,533 & 0,912 & 0,498 & 0,885\\ \mathbf{Pieza\, de\, Tela} & 6 & 7 & 8 & 9 & 10\\ \mathbf{Reflectancia} & 0,291 & 0,805 & 0,436 & 0,868 & 0,712\\ \hline \end{array} \]

Para aplicar la prueba de bondad de ajuste Kolmogorov-Smirnov se debe construir la tabla con la distribución de frecuencias observadas y esperadas, y a partir de estas se obtienen los estadísticos de Kolmogorov-Smirnov \(D_{n}^{+}\), \(D_{n}^{-}\) y \(D_n\), definidos en las ecuaciones (1.54), (1.55) y (1.56), respectivamente.

reflectancia <- c(
  0.608, 0.533, 0.912, 0.498, 0.885, 0.291, 0.805,
  0.436, 0.868, 0.712
)
Xi <- sort(unique(reflectancia))
i <- numeric(length(Xi))
for (j in 1:length(Xi)) {
  i[j] <- sum(reflectancia <= Xi[j])
}
datos7 <- data.frame(
  i = i,
  Xi = Xi
) %>%
  mutate(
    Fn = i / length(reflectancia),
    F0 = punif(q = Xi),
    D_mas = Fn - F0,
    D_menos = F0 - (i - 1) / length(reflectancia)
  )
kableExtra::kbl(
  datos7,
  digits = 3,
  booktabs = TRUE,
  row.names = NA,
  col.names = c(
    "$i$",
    "$x_{\\left( i \\right)}$",
    "$F_{n}^{*}\\left( x_{\\left( i \\right)} \\right)$",
    "$F_{0}\\left( x_{\\left( i \\right)} \\right)$",
    "$D_{n}^{+}$",
    "$D_{n}^{-}$"
  ),
  align = c("c", "c", "c", "c", "c", "c"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  caption = "\\label{tab2:datos7}Distribución empírica de la reflectancia Vs la distribución $U\\left(a=0, b =1 \\right)$",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  column_spec(
    5,
    color = if_else(datos7$D_mas == max(datos7$D_mas), "red", "black"),
    bold = if_else(datos7$D_mas == max(datos7$D_mas), TRUE, FALSE)
  ) %>%
  column_spec(
    6,
    color = if_else(datos7$D_menos == max(datos7$D_menos), "red", "black"),
    bold = if_else(datos7$D_menos == max(datos7$D_menos), TRUE, FALSE)
  ) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 1.17: Distribución empírica de la reflectancia Vs la distribución \(U\left(a=0, b =1 \right)\)

\(i\)	\(x_{\left( i \right)}\)	\(F_{n}^{*}\left( x_{\left( i \right)} \right)\)	\(F_{0}\left( x_{\left( i \right)} \right)\)	\(D_{n}^{+}\)	\(D_{n}^{-}\)
1	0,291	0,1	0,291	-0,191	0,291
2	0,436	0,2	0,436	-0,236	0,336
3	0,498	0,3	0,498	-0,198	0,298
4	0,533	0,4	0,533	-0,133	0,233
5	0,608	0,5	0,608	-0,108	0,208
6	0,712	0,6	0,712	-0,112	0,212
7	0,805	0,7	0,805	-0,105	0,205
8	0,868	0,8	0,868	-0,068	0,168
9	0,885	0,9	0,885	0,015	0,085
10	0,912	1,0	0,912	0,088	0,012

Como se observa en la tabla 1.17 \(D_{n}^{+} = 0,088\) y \(D_{n}^{-} = 0,336\). En consecuencia, el estadístico bilateral de Kolgomoro-Smirnov es \(D_{n} = 0,336\). Mientras que \(p_{valor} = 0,165092\).

Dado que el \(p_{valor}>0,05\) se acepta la hipótesis de que la reflectancia tiene distribución \(U(a=0, b=1)\).

Los resultados anteriores se pueden obtener con el siguiente código.

D_mas <- max(datos7$D_mas)
D_menos <- max(datos7$D_menos)
Dn <- max(D_mas, D_menos)
p.valor <- 1 - pkolmim(d = Dn, n = length(reflectancia))
paste("El valor de estadístico Dn_mas =", D_mas)
paste("El valor de estadístico Dn_menos =", D_menos)
paste("El valor de estadístico Dn =", Dn)
paste("El p.valor =", p.valor)

#> [1] "El valor de estadístico Dn_mas = 0.088"
#> [1] "El valor de estadístico Dn_menos = 0.336"
#> [1] "El valor de estadístico Dn = 0.336"
#> [1] "El p.valor = 0.165092038878279"

La figura 1.2, generada por el código que se presenta luego, compara la gráfica de la función de distribución acumulada empíricas de la muestra y la distribución \(U(a=0, b=1)\).

ggplot(as.data.frame(reflectancia), aes(x = reflectancia)) +
  xlim(0, 1) +
  stat_ecdf(
    aes(colour = "Distribución empírica"),
    geom = "step", color = "steelblue"
  ) +
  stat_ecdf(geom = "point", color = "black") +
  stat_function(
    aes(colour = "Distribución uniforme"),
    fun = punif, color = "red",
    args = list(min = 0, max = 1)
  ) +
  labs(
    x = "Reflectancia",
    y = "Distribución acumulada"
  ) +
  theme(legend.position = "top") +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(
      colour = "lightsteelblue1",
      size = 0.25, linetype = "twodash"
    ),
    axis.title = element_text(face = "bold.italic"),
    panel.background = element_rect(fill = "white", linetype = "solid"),
    plot.background = element_rect(fill = "antiquewhite", colour = NA)
  ) +
  theme(
    axis.line = element_line(
      colour = NA,
      linetype = "solid"
    ), axis.title = element_text(size = 10),
    axis.text.x = element_text(size = 9),
    axis.text.y = element_text(size = 9),
    plot.title = element_text(size = 12),
    plot.background = element_rect(fill = "antiquewhite1")
  )

Figura 1.2: Función de distribución acumulada empírica de la reflectancia

Luego, la ejecución del siguiente código muestra la prueba de bondad de ajuste Kolmogorov-Smirnov, para determinar si la muestra obtenida de la variable reflectancia se puede modelar por la distribución uniforme estándar.

ks.test(reflectancia, "punif",
  min = 0, max = 1,
  alternative = "two.sided"
)

#> 
#>  One-sample Kolmogorov-Smirnov test
#> 
#> data:  reflectancia
#> D = 0.336, p-value = 0.1651
#> alternative hypothesis: two-sided

Note que los valores \(D_n\) y el \(p_{valor}\) obtenidos en esta salida coinciden con los anteriores.

Ejemplo 1.9 (Prueba de Bondad de Ajuste Kolmogorov-Smirnov) La prueba de bondad de ajuste Chi-Cuadrada aplicada en el ejemplo 1.6 mostró que la longitud de la pieza no se distribuyen normalmente con media 10,5 y desviación estándar 0,15. A continuación se procede a verificar esta conclusión aplicando la prueba de bondad de ajuste Kolmogorov-Smirnov.

La tabla 1.18 muestra la función de distribución acumulada empíricas y la función de distribución acumulada propuesta en \(H_0\).

lp <- c(
  10.39, 10.66, 10.12, 10.32, 10.25, 10.91, 10.52, 10.83,
  10.72, 10.28, 10.35, 10.46, 10.54, 10.72, 10.23, 10.18,
  10.62, 10.49, 10.32, 10.61, 10.64, 10.23, 10.29, 10.78,
  10.81, 10.39, 10.34, 10.62, 10.75, 10.34, 10.41, 10.81,
  10.64, 10.53, 10.31, 10.46, 10.47, 10.43, 10.57, 10.74
)
Xi <- sort(unique(lp))
i <- numeric(length(Xi))
for (j in 1:length(Xi)) {
  i[j] <- sum(lp <= Xi[j])
}
datos8 <- data.frame(
  i = i,
  Xi = Xi
) %>%
  mutate(
    Fn = i / length(lp),
    F0 = pnorm(q = Xi, mean = 10.5, sd = 0.15),
    D_mas = Fn - F0,
    D_menos = F0 - (i - 1) / length(lp)
  )
kableExtra::kbl(
  datos8,
  digits = 4,
  booktabs = TRUE,
  row.names = NA,
  col.names = c(
    "$i$",
    "$x_{\\left( i \\right)}$",
    "$F_{n}^{*}\\left( x_{\\left( i \\right)} \\right)$",
    "$F_{0}\\left( x_{\\left( i \\right)} \\right)$",
    "$D_{n}^{+}$",
    "$D_{n}^{-}$"
  ),
  align = c("c", "c", "c", "c", "c", "c"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  caption = "\\label{tab2:datos8}Distribución empírica de la longitud de la pieza Vs la distribución $N\\left(\\mu=10,5;\\, \\sigma =0,15 \\right)$",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  column_spec(
    5,
    color = if_else(datos8$D_mas == max(datos8$D_mas),
                    "red", "black"),
    bold = 
      if_else(datos8$D_mas == max(datos8$D_mas), 
              TRUE, FALSE)
  ) %>%
  column_spec(
    6,
    color = if_else(datos8$D_menos == max(datos8$D_menos), "red", "black"),
    bold = if_else(datos8$D_menos == max(datos8$D_menos), TRUE, FALSE)
  ) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 1.18: Distribución empírica de la longitud de la pieza Vs la distribución \(N\left(\mu=10,5;\, \sigma =0,15 \right)\)

\(i\)	\(x_{\left( i \right)}\)	\(F_{n}^{*}\left( x_{\left( i \right)} \right)\)	\(F_{0}\left( x_{\left( i \right)} \right)\)	\(D_{n}^{+}\)	\(D_{n}^{-}\)
1	10,12	0,025	0,0056	0,0194	0,0056
2	10,18	0,050	0,0164	0,0336	-0,0086
4	10,23	0,100	0,0359	0,0641	-0,0391
5	10,25	0,125	0,0478	0,0772	-0,0522
6	10,28	0,150	0,0712	0,0788	-0,0538
7	10,29	0,175	0,0808	0,0942	-0,0692
8	10,31	0,200	0,1026	0,0974	-0,0724
10	10,32	0,250	0,1151	0,1349	-0,1099
12	10,34	0,300	0,1431	0,1569	-0,1319
13	10,35	0,325	0,1587	0,1663	-0,1413
15	10,39	0,375	0,2317	0,1433	-0,1183
16	10,41	0,400	0,2743	0,1257	-0,1007
17	10,43	0,425	0,3204	0,1046	-0,0796
19	10,46	0,475	0,3949	0,0801	-0,0551
20	10,47	0,500	0,4207	0,0793	-0,0543
21	10,49	0,525	0,4734	0,0516	-0,0266
22	10,52	0,550	0,5530	-0,0030	0,0280
23	10,53	0,575	0,5793	-0,0043	0,0293
24	10,54	0,600	0,6051	-0,0051	0,0301
25	10,57	0,625	0,6796	-0,0546	0,0796
26	10,61	0,650	0,7683	-0,1183	0,1433
28	10,62	0,700	0,7881	-0,0881	0,1131
30	10,64	0,750	0,8247	-0,0747	0,0997
31	10,66	0,775	0,8569	-0,0819	0,1069
33	10,72	0,825	0,9288	-0,1038	0,1288
34	10,74	0,850	0,9452	-0,0952	0,1202
35	10,75	0,875	0,9522	-0,0772	0,1022
36	10,78	0,900	0,9690	-0,0690	0,0940
38	10,81	0,950	0,9806	-0,0306	0,0556
39	10,83	0,975	0,9861	-0,0111	0,0361
40	10,91	1,000	0,9969	0,0031	0,0219

Como se observa en la tabla 1.18 \(D_{n}^{+} = 0,1663447\) y \(D_{n}^{-} = 0,1433224\). En consecuencia, el estadístico bilateral de Kolgomoro-Smirnov es \(D_{n} = 0,1663447\). Mientras que \(p_{valor} = 0,9038855\).

Dado que el \(p_{valor} > 0,05\) se acepta la hipótesis de que la longitud de la pieza tiene distribución \(N\left(\mu=10,5;\; \sigma =0,15 \right)\).

Los resultados anteriores se pueden obtener con el siguiente código.

D_mas <- max(datos8$D_mas)
D_menos <- max(datos8$D_menos)
Dn <- max(D_mas, D_menos)
p.valor <- 1 - pkolmim(d = Dn, n = length(lp))
paste("El valor de estadístico Dn_mas =", D_mas)
paste("El valor de estadístico Dn_menos =", D_menos)
paste("El valor de estadístico Dn =", Dn)
paste("El p.valor =", p.valor)

#> [1] "El valor de estadístico Dn_mas = 0.166344746068544"
#> [1] "El valor de estadístico Dn_menos = 0.143322425365201"
#> [1] "El valor de estadístico Dn = 0.166344746068544"
#> [1] "El p.valor = 0.195116738373406"

La figura 1.3, generada por el código que se presenta luego, compara la gráfica de la función de distribución acumulada empíricas y la distribución \(N\left(\mu=10,5;\; \sigma =0,15 \right)\). Como se nota en esta gráfica, la distribución empírica se ajusta a la distribución \(N\left(\mu=10,5;\; \sigma =0,15 \right)\), representada por la curva roja.

ggplot(as.data.frame(lp), aes(x = lp)) +
  stat_ecdf(
    aes(colour = "Distribución empírica"),
    geom = "step", color = "steelblue"
  ) +
  stat_ecdf(geom = "point", color = "black") +
  stat_function(
    aes(colour = "Distribución Normal"),
    fun = pnorm, color = "red",
    args = list(mean = 10.5, sd = 0.15)
  ) +
  labs(
    x = "Longitud de la pieza",
    y = "Distribución acumulada"
  ) +
  theme(legend.position = "top") +
  theme_minimal() +
  theme(
    panel.grid.major = element_line(
      colour = "lightsteelblue1",
      size = 0.25, linetype = "twodash"
    ),
    axis.title = element_text(face = "bold.italic"),
    panel.background = 
      element_rect(fill = "white", linetype = "solid"),
    plot.background = 
      element_rect(fill = "antiquewhite", colour = NA)
  ) +
  theme(
    axis.line = element_line(
      colour = NA,
      linetype = "solid"
    ), axis.title = element_text(size = 10),
    axis.text.x = element_text(size = 9),
    axis.text.y = element_text(size = 9),
    plot.title = element_text(size = 12),
    plot.background = element_rect(fill = "antiquewhite1")
  )

Figura 1.3: Función de distribución acumulada empírica de la longitud de la pieza

Luego, la ejecución del siguiente código muestra la prueba de bondad de ajuste Kolmogorov-Smirnov, para determinar si la muestra obtenida de la longitud de la pieza se puede modelar por la distribución \(N\left(\mu=10,5;\; \sigma =0,15 \right)\).

ks.test(lp, "pnorm",
  mean = 10.5, sd = 0.15,
  alternative = "two.sided", exact = TRUE
)

#> 
#>  One-sample Kolmogorov-Smirnov test
#> 
#> data:  lp
#> D = 0.16634, p-value = 0.1951
#> alternative hypothesis: two-sided

Note que los valores \(D_n\) y \(p_{valor}\) son iguales a los obtenidos anteriormente.

1.5.3 Problemas Propuestos

1. Un estudiante de estadística desea verificar si es razonable suponer que algunos datos sobre ventas han sido tomados de una población normal antes de llevar a cabo una prueba de hipótesis sobre la media de las ventas. Recoge algunos datos sobre ventas, calcula \(\bar{x}= 78\) y \(s = 9\), y tabula los datos de la siguiente manera: \[ \begin{array}{lcccccc} \hline \mathbf{Nivel \, de\, ventas} & \leq65 & 66 - 70 & 71 - 75 & 76 - 80 & 81 - 85 & \geq 86 \\ \mathbf{N^{\circ} \, de\, observaciones} & 10 & 20 & 40 & 50 & 40 & 40\\ \hline \end{array} \]
   1. ¿Es importante para el estudiante de estadística verificar si los datos están distribuidos normalmente? Explique su respuesta.
   2. Al nivel de significancia de 5%. ¿la distribución de frecuencias observada sigue una distribución normal \(\bar{x}= 78\) y \(s = 9\)?
2. Se utiliza un proceso de plateado para cubrir cierto tipo de charola de servicio. Cuando el proceso está bajo control el espesor de la plata sobre la charola variara de forma aleatoria siguiendo una distribución normal con una media de 0,02 milímetros y una desviación estándar de 0,005 milímetros. Suponga que las siguientes 12 charolas examinadas muestran los siguientes espesores de plata: 0,019; 0,021; 0,020; 0,019; 0,020; 0,018; 0,023; 0,021; 0,024; 0,022; 0,023; 0,022. Utilice la prueba de Kolmogorov-Smirnov para determinar si está bajo control. Utilice \(\alpha= 0,05\).

1.6 Prueba de Rachas o Carreras de Walf-Wolfowits

Si un investigador desea llegar a alguna conclusión acerca de una población usando la información contenida en una muestra extraída de esa población, entonces la muestra debe ser aleatoria; es decir, las observaciones sucesivas deben ser independientes. En este sentido la prueba de la racha o carrera se usa para probar la hipótesis de que una muestra tomada de una población es aleatoria.

1.6.1 Racionalización de la Prueba de Rachas o Carreras de Walf-Wolfowits

La prueba de la racha está basada en el número de rachas, corridas, carreras o series que exhibe una muestra. Una racha se define como una sucesión de símbolos idénticos que son seguidos y precedidos por diferentes símbolos o por ningún símbolo. Por ejemplo, supóngase que se tiene una población que se ha codificada en forma binaria, es decir solo puede tomar los valores 0 ó 1, y que se presentan las muestras

\[ \begin{array}{cc} 00000000111111 & \acute{o} & 10101010101010 \\ \end{array} \]

no cabe duda de que tienen un aspecto poco aleatorio y que cabe dudar de que se hayan observado en condiciones aleatorias. En cambio, la sospecha no surgiría, por ejemplo, ante las muestras:

\[ \begin{array}{cc} 10011100110001 & \acute{o} & 01110010000111 \\ \end{array} \]

Ante el desconocimiento de la distribución poblacional (en este caso, de las proporciones de 0 y 1 que debería presentarse) el único rasgo por el que se puede diferenciar la poca aleatoriedad de las primeras respecto de las segundas es por la colocación de los dígitos y, más exactamente, por el número, \(R\), de rachas de cada uno que haya presentes. Así, en la primera es \(R=2\), en la segunda \(R=14\), es \(R=7\) y \(R=6\) en las dos últimas.

Naturalmente son más las secuencias con un número intermedio de rachas que con un número extremo, de forma que una secuencia elegida al azar contendrá, muy probablemente, un número intermedio de rachas. En tal sentido, regiones críticas de la forma

\[ \left\{ R< k_{1}\right\} \cup \left\{R> k_{2} \right\} \]

sirven entonces como test para el contraste de la hipótesis de que la muestra es aleatoria simple.

Suponga de manera general, que en la muestra de tamaño \(n\) hay \(n_1\) elementos de una clase y \(n_2\) elementos de otra clase, es decir la muestra consta de \(n =n_1 + n_2\) elementos binarios. Entonces, la distribución de \(R\), bajo \(H_0\), se puede hallar explícitamente por:

\[ \begin{equation} \left.\begin{matrix} P\left\{ R=2r\right\}=2 \frac{\binom{n_{1} - 1}{r-1} \binom{n_{2}-1}{r-1}}{\binom{n}{n_{1}}} \\ P\left\{ R=2r +1\right\}= \frac{\binom{n_{1} - 1}{r-1} \binom{n_{2}-1}{r}+\binom{n_{1}-1}{r}\binom{n_{2}-1}{r-1}}{\binom{n}{n_{1}}} \end{matrix}\right\}para \: r \: mín\left ( n_{1},n_{2} \right ). \tag{1.57} \end{equation} \]

En efecto, \(n_1\) elementos de una clase y \(n_2\) elementos de la otra clase pueden ordenarse de \((n_1+n_2 )!=n!\) maneras, igualmente probables. Para formar secuencias con \(r\) rachas de elementos de una clase y \(r\) rachas de la otra clase, los elementos de la primera clase se pueden ordenar de \(n_1!\) formas, dividiéndose luego \(r\) en grupos, mediante la elección de \(r-1\) de los \(n_1-1\) huecos existentes entre ellos.

Análogamente, los elementos de la segunda clase pueden ser ordenados de \(n_2 !\) maneras y divididos en \(r\) grupos de \(\binom{n_{2}-1}{r-1}\) formas. Después hay que intercambiar los grupos formados, empezando con el primer grupo de elementos de la primera clase, o bien con el primer grupo de elementos de la segunda clase. En total, la probabilidad de que haya \(r\) rachas de cada tipo es

\[ \frac{n!\binom{n_1 -1}{r-1}n_2!\binom{n_2-1}{r-1}}{\left ( n_1 + n_2 \right )!}=2 \frac{\binom{n_{1} - 1}{r-1} \binom{n_{2}-1}{r-1}}{\binom{n}{n_{1}}}. \]

De manera similar se calcula la probabilidad de que haya \(r\) rachas de elementos de la primera clase y \(r+1\) rachas de elementos de la segunda clase y la probabilidad de que haya \(r+1\) rachas de elementos de la primera clase y \(r\) rachas de elementos de la segunda clase, las cuales dan lugar a los dos sumandos de la segunda ecuación de la expresión (1.57). (Naturalmente, el número de rachas de un tipo y otro se diferencian, a lo sumo, en 1).

Si \(n_1\) y \(n_2\) son grandes (superiores a 20), puede probarse que

\[ R\simeq N\left ( \mu_R =\frac{2n_1n_2}{n} +1,\sigma_R=\sqrt{\frac{2n_1n_2\left ( 2n_1n_2-n \right )}{n^{2}\left ( n-1 \right )}} \right ). \]

De donde se establece que el estadístico

\[ \begin{equation} Z=\frac{R+h-\frac{2n_1n_2}{n} -1}{\sqrt{\frac{2n_1n_2\left ( 2n_1n_2-n \right )}{n^{2}\left ( n-1 \right )}}}\simeq N\left(\mu_Z=0, \sigma_Z=1 \right). \tag{1.58} \end{equation} \]

Donde \(h=+0,5\) si \(r<(2n_1 n_2)⁄n \;+1\), y \(h=-0,5\) si \(r>(2n_1 n_2)⁄n\;+1\).

Para contrastar la hipótesis nula \(H_0:\) la cual representa que la muestra es aleatoria usando el estadístico \(R\) es necesario determinar el valor crítico de \(R\) para un nivel de significancia \(\alpha\) especificado y los valores de los parámetros \(n_1\) y \(n_2\), determinados a partir de la muestra. En la literatura existen diversas tablas que proporcionan estos valores, generalmente para \(n_1, n_2 \leq20\). Otra forma de contrastar \(H_0\) es a través del \(p_{valor}\) calculado por medio de la función de distribución acumulada de \(R\).

El paquete randtests proporciona las siguientes funciones relacionadas con la distribución de las rachas:

druns(x, n1, n2, log = FALSE)
pruns(q, n1, n2, lower.tail = TRUE, log.p = FALSE)
qruns(p, n1, n2, lower.tail = TRUE, log.p = FALSE)
rruns(n, n1, n2)

La función druns determina la función de masa de probabilidad de \(R\), la función pruns evalúa la función de distribución de \(R\), la función qruns determina valores críticos de \(R\) y rruns genera números aleatorios de \(R\). Los parámetros de estas funciones se describen en la tabla ??.

knitr::kable(
  data.frame(
    stringsAsFactors = FALSE,
    Parámetros = c(
      "x, q", "p", "n",
      "n1, n2", "log, log.p", "lower.tail"
    ),
    Descripción = c(
      "Vector numérico de cuantiles", "Vector numérico de probabilidades",
      "Número valores aleatorios de rachas a generar",
      "El número de elementos del primer y segundo tipo, respectivamente",
      "Lógico; si es TRUE, las probabilidades p se dan como log(p)",
      "Lógico; si es TRUE (por defecto), las probabilidades son $P[X \\leq x]$, de lo contrario, $P[X > x]$."
    )
  ),
  booktabs = TRUE,
  caption = "Descripción de los parámetros de funciones de R relacionadas con la distribución de las rachas",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2()

Tabla 1.19: Descripción de los parámetros de funciones de R relacionadas con la distribución de las rachas

Parámetros	Descripción
x, q	Vector numérico de cuantiles
p	Vector numérico de probabilidades
n	Número valores aleatorios de rachas a generar
n1, n2	El número de elementos del primer y segundo tipo, respectivamente
log, log.p	Lógico; si es TRUE, las probabilidades p se dan como log(p)
lower.tail	Lógico; si es TRUE (por defecto), las probabilidades son \(P[X \leq x]\), de lo contrario, \(P[X > x]\).

Por otro lado, la función runs.test(x, alternative, threshold, pvalue, plot) del paquete randtests, permite aplicar la prueba de las rachas de Walf-Wolfowits. Los parámetros de esta función se especifican en la tabla 1.20.

knitr::kable(
  data.frame(
    stringsAsFactors = FALSE,
    Parámetros = c("x", "alternative", "threshold", "pvalue", "plot"),
    Descripción = c(
      "Vector numérico que contiene las observaciones",
      "Cadena de caracteres indicando la hipótesis alternativa. Puede tomar los valores: \"two.sided\" (por defecto), \"left.sided\" o \"right.sided\". Puede especificar sólo la letra inicial",
      "El punto de corte para transformar los datos en un vector dicotómico. Por defecto toma el valor de la mediana de x",
      "Cadena de caracteres que especifica el método utilizado para calcular el valor p. Puede ser normal (por defecto), o exacto",
      "Valor lógico para indicar si se debe crear un gráfico. Si es 'TRUE', entonces el gráfico será graficado"
    )
  ),
  booktabs = TRUE,
  caption = "\\label{tab2:par-runstest}Descripción de los parámetros de función `runs.test`",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2()

Tabla 1.20: Descripción de los parámetros de función runs.test

Parámetros	Descripción
x	Vector numérico que contiene las observaciones
alternative	Cadena de caracteres indicando la hipótesis alternativa. Puede tomar los valores: “two.sided” (por defecto), “left.sided” o “right.sided”. Puede especificar sólo la letra inicial
threshold	El punto de corte para transformar los datos en un vector dicotómico. Por defecto toma el valor de la mediana de x
pvalue	Cadena de caracteres que especifica el método utilizado para calcular el valor p. Puede ser normal (por defecto), o exacto
plot	Valor lógico para indicar si se debe crear un gráfico. Si es ‘TRUE’, entonces el gráfico será graficado

1.6.2 Implementación en R de la Prueba de Rachas o Carreras de Walf-Wolfowits

Ejemplo 1.10 (Prueba de Rachas o Carreras de Walf-Wolfowits) En un estudio de la dinámica de agresión en niños pequeños, un investigador observó pares de niños en una situación de juego controlada. La mayoría de los 24 niños que sirvieron como sujetos en el estudio provenía de la misma guardería y, por lo tanto, jugaban juntos diariamente. Ya que el observador fue capaz de ingeniarse para observar sólo dos niños en cualquier día, estaba interesado en que podría haberse introducido sesgo en el estudio por discusiones entre aquellos niños que ya habían servido como sujetos y aquellos que sirvieron posteriormente. Si tales discusiones tenían algún efecto en el nivel de agresión en las sesiones de juego, este efecto podría mostrarse como una carencia de aleatoriedad en las puntuaciones de agresión en el orden en que fueron colectadas. Después de concluir el estudio, la aleatoriedad de la secuencia de puntuaciones fue probada al convertir las puntuaciones de agresión de cada niño a un signo más o a un signo menos, dependiendo si se encontraba por arriba o por debajo de la mediana del grupo. Los datos se muestran en la tabla 1.21.

puntuacion <- c(
  31L, 23L, 36L, 43L, 51L, 44L, 12L, 26L, 43L, 75L, 2L, 3L,
  15L, 18L, 78L, 24L, 13L, 27L, 86L, 61L, 13L, 7L, 6L, 8L
)
datos9 <- data.frame(
  niños = 1:length(puntuacion),
  puntuacion = puntuacion
) %>%
  mutate(
    signo = as.factor(
      case_when(
        puntuacion > median(puntuacion, na.rm = TRUE) ~ "\\+",
        puntuacion == median(puntuacion, na.rm = TRUE) ~ "=",
        puntuacion < median(puntuacion, na.rm = TRUE) ~ "\\-"
      )
    )
  )
kableExtra::kbl(
  datos9,
  digits = 4,
  booktabs = TRUE,
  row.names = NA,
  col.names = c(
    "Niños",
    "Puntuación",
    "Posición de la puntuación con respecto a la mediana "
  ),
  align = c("c", "c", "c"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  caption = "\\label{tab2:datos9}Puntuaciones de agresión de acuerdo al
orden de ocurrencia",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 1.21: Puntuaciones de agresión de acuerdo al orden de ocurrencia

Niños	Puntuación	Posición de la puntuación con respecto a la mediana
1	31	+
2	23	-
3	36	+
4	43	+
5	51	+
6	44	+
7	12	-
8	26	+
9	43	+
10	75	+
11	2	-
12	3	-
13	15	-
14	18	-
15	78	+
16	24	-
17	13	-
18	27	+
19	86	+
20	61	+
21	13	-
22	7	-
23	6	-
24	8	-

Los datos fueron tomados de Siegel (2012, pág. 85).

El siguiente trozo de código permite determinar el valor del estadístico \(R\), el \(p_{valor}\), la media y la varianza de \(R\).

mediana <- datos9 %>%
  pull(puntuacion) %>%
  median(., na.rm = TRUE)
n1 <- datos9 %>%
  dplyr::filter(signo == "\\+") %>%
  pull(signo) %>%
  as.character(.) %>%
  length(.)
n2 <- datos9 %>%
  dplyr::filter(signo == "\\-") %>%
  pull(signo) %>%
  as.character(.) %>%
  length(.)
n <- n1 + n2
rachas <- rle(sign(na.omit(puntuacion) - mediana))
rachas9 <- as.data.frame.list(rachas) %>%
  dplyr::rename(
    lonjitud = lengths,
    signo = values
  )
r1 <- rachas9 %>%
  dplyr::filter(
    signo == 1
  ) %>%
  nrow(.)
r2 <- rachas9 %>%
  dplyr::filter(
    signo == -1
  ) %>%
  nrow(.)
r <- r1 + r2
mu_R <- 2 * n1 * n2 / n + 1
var_R <- 2 * n1 * n2 * (2 * n1 * n2 - n) / (n^2 * (n - 1))
p.valor <- 2 * pruns(q = r, n1 = n1, n2 = n2)
paste("La media del estadístico R =", mu_R)
paste("La varianza del estadístico R =", var_R)
paste("El valor del estadístico R =", r)
paste("El p.valor =", round(x = p.valor, digits = 4))

#> [1] "La media del estadístico R = 13"
#> [1] "La varianza del estadístico R = 5.73913043478261"
#> [1] "El valor del estadístico R = 10"
#> [1] "El p.valor = 0.3009"

Note en la salida anterior que el valor del estadístico \(R = 10\). Mientras que, el \(p_{valor} = 0,300889\), por lo tanto las puntuaciones de agresión de acuerdo al orden en que fueron colectadas tienen un comportamiento aleatorio, lo que implica que las discusiones no tienen efecto en el nivel de agresión en las sesiones de juego.

La figura 1.4 muestra la gráfica de las rachas por debajo y por arriba de la mediana de la muestra (linea horizontal roja). Y en esta se observa que las longitudes de las rachas tienen un comportamiento aleatorio, lo cual contrasta con la conclusión anterior.

ggplot(datos9, aes(x = niños, y = puntuacion)) +
  geom_point(aes(colour = as.factor(signo)), size = 2) +
  geom_vline(
    xintercept = rachas9 %>%
      dplyr::select(lonjitud) %>%
      dplyr::slice(1:nrow(.) - 1) %>%
      dplyr::pull(lonjitud) %>%
      cumsum(.) + 0.5,
    color = "steelblue",
    linetype = "dashed",
    size = 1
  ) +
  geom_hline(
    yintercept = mediana,
    color = "red",
    linetype = "solid",
    size = 1
  ) +
  labs(
    x = "Niños",
    y = "Puntuación"
  ) +
  theme_minimal() +
  scale_colour_brewer(palette = "Set1") +
  theme(
    panel.grid.major = element_line(
      colour = "lightsteelblue1",
      size = 0.25, linetype = "twodash"
    ),
    axis.title = element_text(face = "bold.italic"),
    panel.background = 
      element_rect(fill = "white", linetype = "solid"),
    plot.background = 
      element_rect(fill = "antiquewhite", colour = NA)
  ) +
  theme(
    axis.line = element_line(
      colour = NA,
      linetype = "solid"
    ), axis.title = element_text(size = 10),
    axis.text.x = element_text(size = 9),
    axis.text.y = element_text(size = 9),
    plot.title = element_text(size = 12),
    plot.background = element_rect(fill = "antiquewhite1"),
    legend.position = "NULL"
  )

Figura 1.4: Rachas para las puntuaciones de agresión de acuerdo al orden de ocurrencia

R tiene implementadas muchas funciones que ejecutan la prueba de Rachas o Carreras de Walf-Wolfowits. Una de ellas es la función runs.test del paquete randtests, cuyos parámetros se definen en la tabla 1.20. El siguiente script ejecuta la prueba para contrastar la hipótesis de aleatoriedad de la muestra.

randtests::runs.test(
  x = puntuacion, alternative = "two.sided",
  threshold = median(puntuacion), pvalue = "exact",
  plot = FALSE
)

#> 
#>  Runs Test
#> 
#> data:  puntuacion
#> statistic = -1.2523, runs = 10, n1 = 12, n2 = 12, n = 24, p-value =
#> 0.3009
#> alternative hypothesis: nonrandomness

La función runs.test del paquete randtests también realiza el gráfico mostrado en la figura 1.4, pero menos elegante que este. Si se quiere ver este gráfico se cambia el parámetro plot a TRUE.

La función runs.test del paquete snpar también ejecuta la prueba, como se muestra en el siguiente trozo de código.

snpar::runs.test(
  x = puntuacion, exact = T,
  alternative = "two.sided"
)

#> 
#>  Exact runs test
#> 
#> data:  puntuacion
#> Runs = 10, p-value = 0.3009
#> alternative hypothesis: two.sided

Ejemplo 1.11 (Prueba de Rachas o Carreras de Walf-Wolfowits) Un investigador estaba interesado en averiguar si la disposición de hombres y mujeres en una fila enfrente de la taquilla de un teatro era un arreglo aleatorio. Los datos se obtuvieron simplemente anotando el sexo de cada una de las 50 personas al aproximarse a la taquilla. \[ \begin{array}{ccccccccccccc} \hline M & F & M & F & M & M & M & F & F & M & F & M & F \\ M & F & M & M & M & M & F & M & F & M & F & M & M \\ F & F & F & M & F & M & F & M & F & M & M & F & \\ M & M & F & M & M & M & M & F & M & F & M & M & \\ \hline \end{array} \] Los datos fueron tomados de Siegel (2012, pág. 86)

cola <- c(
  "M", "F", "M", "F", "M", "M", "M", "F", "F", "M", "F",
  "M", "F", "M", "F", "M", "M", "M", "M", "F", "M", "F",
  "M", "F", "M", "M", "F", "F", "F", "M", "F", "M", "F",
  "M", "F", "M", "M", "F", "M", "M", "F", "M", "M", "M",
  "M", "F", "M", "F", "M", "M"
)
n1 <- min(table(cola))
n2 <- max(table(cola))
n <- n1 + n2
ind <- fct_infreq(cola) %>%
  fct_rev() %>%
  levels(.)
n1 <- length(cola[cola == ind[1]])
n2 <- length(cola[cola == ind[2]])
rachas <- rle(cola)
rachas10 <- as.data.frame.list(rachas) %>%
  dplyr::rename(
    lonjitud = lengths,
    signo = values
  )
r1 <- rachas10 %>%
  dplyr::filter(
    signo == ind[1]
  ) %>%
  nrow(.)
r2 <- rachas10 %>%
  dplyr::filter(
    signo == ind[2]
  ) %>%
  nrow(.)
r <- r1 + r2
mu_R <- 2 * n1 * n2 / n + 1
var_R <- 2 * n1 * n2 * (2 * n1 * n2 - n) / (n^2 * (n - 1))
z <- (r - mu_R) / sqrt(var_R)
z.corr <- (r - mu_R - 0.5) / sqrt(var_R)
p.valor <- 2 * pruns(q = r - 1, n1 = n1, n2 = n2, 
                     lower.tail = FALSE)
p.valor.aprox <- 2 * pnorm(z, lower.tail = FALSE)
p.valor.aprox.cor <- 2 * pnorm(z.corr, lower.tail = FALSE)
paste("n =", n, "n1 =", n1, "n2 =", n2, "r1 =", r1, 
      "r2 =", r2)
paste("La media del estadístico R =", mu_R)
paste("La varianza del estadístico R =", var_R)
paste("El valor del estadístico R =", r)
paste("El valor del estadístico Z =", z)
paste("El valor del estadístico Z (corregido) =", z.corr)
paste("El p.valor =", round(x = p.valor, digits = 6))
paste("El p.valor.aprox =", round(x = p.valor.aprox, digits = 6))
paste("El p.valor.aprox con corrección =", round(x = p.valor.aprox.cor, digits = 6))

#> [1] "n = 50 n1 = 20 n2 = 30 r1 = 17 r2 = 18"
#> [1] "La media del estadístico R = 25"
#> [1] "La varianza del estadístico R = 11.265306122449"
#> [1] "El valor del estadístico R = 35"
#> [1] "El valor del estadístico Z = 2.97939785765562"
#> [1] "El valor del estadístico Z (corregido) = 2.83042796477284"
#> [1] "El p.valor = 0.003748"
#> [1] "El p.valor.aprox = 0.002888"
#> [1] "El p.valor.aprox con corrección = 0.004649"

Como se observa en la salida anterior el \(p_{valor} = 0,003748\), lo que indica que la disposición de hombres y mujeres en la fila enfrente de la taquilla del teatro no es aleatorio. Note que \(R = 35 > \mu_R = 25\), por lo que el \(p_{valor}=2P\left(R \geq r \right)\).

La figura 1.5 muestra la gráfica de las rachas de hombres y mujeres en la taquilla. Y esta muestra que las longitudes de las rachas no tienen un comportamiento aleatorio, lo cual contrasta con la conclusión anterior.

datos10 <- tibble(
  personas = 1:length(cola),
  cola = cola,
  codigo = ifelse(cola == "M", 1, 0)
)
ggplot(datos10, aes(x = personas, y = codigo)) +
  geom_point(aes(colour = as.factor(cola)), size = 2) +
  geom_vline(
    xintercept = rachas10 %>%
      dplyr::select(lonjitud) %>%
      dplyr::slice(1:nrow(.) - 1) %>%
      dplyr::pull(lonjitud) %>%
      cumsum(.) + 0.5,
    color = "steelblue",
    linetype = "dashed",
    size = 1
  ) +
  geom_hline(
    yintercept = 0.5,
    color = "red",
    linetype = "solid",
    size = 1
  ) +
  labs(
    x = "Niños",
    y = "Puntuación"
  ) +
  theme_minimal() +
  scale_colour_brewer(palette = "Set1") +
  theme(
    panel.grid.major = element_line(
      colour = "lightsteelblue1",
      size = 0.25, linetype = "twodash"
    ),
    axis.title = element_text(face = "bold.italic"),
    panel.background = element_rect(fill = "white", linetype = "solid"),
    plot.background = element_rect(fill = "antiquewhite", colour = NA)
  ) +
  theme(
    axis.line = element_line(
      colour = NA,
      linetype = "solid"
    ), axis.title = element_text(size = 10),
    axis.text.x = element_text(size = 9),
    axis.text.y = element_text(size = 9),
    plot.title = element_text(size = 12),
    plot.background = element_rect(fill = "antiquewhite1"),
    legend.position = "NULL"
  )

Figura 1.5: Rachas para las puntuaciones de agresión de acuerdo al orden de ocurrencia

Los resultados anteriores se pueden contrastar aplicando la función runs.test del paquete randtests. Note que en el trozo de código que se muestra abajo el parámetro pvalue se ha establecido como "exact" para generar el \(p_{valor}\) exacto. Si este se establece como "normal", el \(p_{valor}\) que muestra la salida es el aproximado sin corrección por continuidad. Intente correr este script cambiando la configuración de este parámetro y observe el \(p_{valor}\).

randtests::runs.test(
  x = datos10$codigo, alternative = "two.sided",
  threshold = 0.5, pvalue = "exact",
  plot = FALSE
)

#> 
#>  Runs Test
#> 
#> data:  datos10$codigo
#> statistic = 2.9794, runs = 35, n1 = 30, n2 = 20, n = 50, p-value =
#> 0.003748
#> alternative hypothesis: nonrandomness

Con el siguiente trozo de código se ha ejecutado la prueba con la función runs.test del paquete snpar. Observe que el parámetro exact = FALSE. Compare este \(p_{valor}\) con los anteriores.

snpar::runs.test(
  x = datos10$codigo, exact = FALSE,
  alternative = "two.sided"
)

#> 
#>  Approximate runs rest
#> 
#> data:  datos10$codigo
#> Runs = 35, p-value = 0.002888
#> alternative hypothesis: two.sided

1.6.3 Problemas Propuestos

1. En una cola hay 20 hombres y 13 mujeres formados para comprar boletos de entrada al juego del Mundialito de Béisbol. Las posiciones fueron: \[ \begin{array}{ccccccccccccc} \hline H & M & H & M & H & H & H & M & H & M & H & H \\ H & M & M & H & H & H & M & M & H & M & H & H \\ H & H & M & H & H & H & M & M & M\\ \hline \end{array} \] ¿Es el proceso aleatorio?

2. Un entrenador de fútbol piensa que el ganar un campeonato un año incrementa la motivación del equipo para ganar el siguiente año. Expresó esta teoría a un estudiante de estadística, quien le pidió los registros de los éxitos y fracasos del equipo de los últimos años. El entrenador le dio una lista, especificando si el equipo había ganado (G) o perdido (P) el campeonato ese año. Los resultados de estas marcas son: \[ \begin{array}{cccccccccc} \hline G & G & G & G & G & G & P & G & G & G \\ G & G & P & G & G & G & G & P & P & G \\ G & G & G & G & G\\ \hline \end{array} \]

1. A un nivel de significancia de 10%. ¿es aleatoria la ocurrencia de éxitos y fracasos?
2. ¿Su respuesta a la parte a), combinada con una inspección ocular de los datos, le dice algo sobre la prueba de corridas?

3. Supóngase que los siguientes datos constituyen la secuencia de residuos para una ecuación de regresión estimada: \[ \begin{array}{ccccccccc} \hline -2,98 & -4,19 & -0,51 & 5,19 & 2,38 & 6,73 & 0,93 & 1,29 & -3,18\\ -1,14 & -0,54 & -2,76 & -1,89 & -4,28 & -0,18 & 0,32 & 0,48 & 1,48\\ -2,43 & -4,69 & 3,18 & 0,64 & 0,89 & 2,08 & 0,98 & -3,28\\ \hline \end{array} \] ¿Existe alguna razón para creer que esta secuencia de residuos no es aleatoria? Úsese \(\alpha=0,05\). Los datos fueron tomados de Canavos (1988, pág. 590).

4. En una línea de producción industrial los artículos se inspeccionan de forma periódica en busca de defectos. La siguiente es una secuencia de artículos defectuosos, D, y no defectuosos, N, producidos por esta línea: \[ \begin{array}{cccccccccccccc} \hline D & D & N & N & N & D & N & N & D & D & N & N & N & N\\ N & D & D & D & N & N & D & N & N & N & N & D & N & D\\ \hline \end{array} \] Utilice la teoría de muestras grandes para la prueba de rachas a un nivel de significancia de 0.05 para determinar si los artículos defectuosos ocurren al azar. Tomado de Walpole (2012, pág. 677).

5. Se utiliza un proceso de plateado para cubrir cierto tipo de charola de servicio. Cuando el proceso está bajo control el espesor de la plata sobre la charola variara de forma aleatoria siguiendo una distribución normal con una media de 0,02 milímetros y una desviación estándar de 0,005 milímetros. Suponga que las siguientes 12 charolas examinadas muestran los siguientes espesores de plata: 0,019; 0,021; 0,020; 0,019; 0,020; 0,018; 0,023; 0,021; 0,024; 0,022; 0,023; 0,022. Utilice la prueba de rachas para determinar si las fluctuaciones en el espesor de una charola a otra son aleatorias. Utilice \(\alpha= 0.05\). Tomado de Walpole (2012, pág. 677).

2 Pruebas que Usan Datos de Dos Muestras Relacionadas

1. Prueba del Signo* (Sección 2.1)

2. Prueba de Rangos Asignados de Wilcoxon (Sección 2.2)

3. Prueba para Datos en Forma de Frecuencias de Mcnemar (Sección 2.3)

2.1 Prueba del Signo

2.1.1 Racionalización de la Prueba del Signo

2.1.2 Implementación en R de la Prueba del Signo

Ejemplo 2.1 (Prueba del Signo) Un investigador estaba estudiando el proceso de toma de decisión esposo-esposa. Se estudió exhaustivamente una muestra de parejas esposo-esposa para determinar el papel percibido de cada uno de ellos respecto de mejorar las adquisiciones domésticas. En cada ocasión, una pareja (cada uno por separado) contestaba un cuestionario concerniente a la influencia que creía ejercer cuando el matrimonio enfrentaba una situación en la que tenía que decidirse la adquisición de enseres para el hogar. Las respuestas a las preguntas se evaluaban mediante una escala que iba de esposo dominante a esposa dominante. Para cada pareja, la diferencia entre sus “percepciones” era determinada y codificada como + si a juicio del esposo, la esposa no debería tener una mayor influencia que él y esto no coincidía con lo informado por la esposa (esposo: “mi opinión debería tener mayor peso que la de ella”, y esposa: “ambos deberíamos ponernos de acuerdo para decidir”). La diferencia se codificaba como (-) cuando ocurría el caso contrario. La diferencia se codificaba como 0 (cero) si la pareja estaba en completo de acuerdo acerca del grado de influencia ejercida en la decisión. Los datos fueron tomados de Siegel (1995, pág. 107). Estos se presentan en la siguiente tabla, por medio del siguiente código.

pareja <- c(LETTERS[1:17])
esposo <- c(
  5L, 4L, 6L, 6L, 3L, 2L, 5L, 3L, 1L, 4L, 5L, 4L, 4L, 7L,
  5L, 5L, 5L
)
esposa <- c(
  3L, 3L, 4L, 5L, 3L, 3L, 2L, 3L, 2L, 3L, 2L, 2L, 5L, 2L,
  5L, 3L, 1L
)
pareja <- tibble(
  pareja = pareja,
  esposo = esposo,
  esposa = esposa,
) %>% 
  mutate(
    signo = as.factor(
      case_when(
        (esposo - esposa) > 0 ~ "\\+",
        (esposo - esposa) == 0 ~ "0",
        (esposo - esposa) < 0 ~ "\\-"
      )
    )
  )
kableExtra::kbl(
  pareja,
  col.names = c("Pareja", "Esposo", "Esposa", "Signo"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  booktabs = TRUE,
  caption = "\\label{tab2:pareja}Juicios acerca de la influencia en la toma de decisiones",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  row_spec(which(pareja$signo == "0"), color = "red",
           strikeout = T) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 2.1: Juicios acerca de la influencia en la toma de decisiones

Pareja	Esposo	Esposa	Signo
A	5	3	+
B	4	3	+
C	6	4	+
D	6	5	+
E	3	3	0
F	2	3	-
G	5	2	+
H	3	3	0
I	1	2	-
J	4	3	+
K	5	2	+
L	4	2	+
M	4	5	-
N	7	2	+
O	5	5	0
P	5	3	+
Q	5	1	+

Aquí la hipótesis nula que se quieres contrastar es que: los esposos y las esposas están de acuerdo en el grado de influencia que ambos deben tener cuando deciden sobre las adquisiciones domésticas. Mientras que la hipótesis alternativa es: los esposos juzgan que ellos deben tener mayor influencia que sus esposas acerca de las decisiones de adquirir enseres para el hogar. Para contrastar estas hipótesis se procede a aplicar la prueba del signo con un nivel de confianza del 95%.

Bajo el supuesto planteado en la hipótesis nula, es de esperarse que el número de diferencias positivas \(S\) se distribuya binomialmente con \(n = 14\) (número de diferencias distintas de cero) y \(p = 1/2\) (probabilidad de que ocurra una diferencia positiva en la muestra), es decir, \(S\sim B(n=14, p=1/2)\). El siguiente código muestra la ejecución de la prueba para los datos previamente mencionados.

n <- pareja %>% 
  dplyr::filter(signo != "0") %>%
  nrow()
s <- pareja %>% 
  dplyr::filter(signo == "\\+") %>%
  nrow()
p.valor <- 1 - pbinom(q = s - 1, size = n, prob = 0.5)
paste("El valor del estadístico S =", s)
paste("El p.valor =", round(x = p.valor, digits = 4))

#> [1] "El valor del estadístico S = 11"
#> [1] "El p.valor = 0.0287"

Como se nota de la salida anterior, el \(p_{valor}=0,028687\), lo que indica que los esposos creen que son ellos los que deben tener una mayor influencia al momento de tomar decisiones acerca de la adquisición de enseres domésticos, en comparación con la que deben tener las esposas. Los resultados anteriores pueden ser obtenidos con la función SIGN.test del paquete BSDA, como se meustra a continuación, con el siguiente segmento de código.

SIGN.test(
  x = esposo, y = esposa, md = 0, alternative = "g",
  conf.level = 0.95
)

#> 
#>  Dependent-samples Sign-Test
#> 
#> data:  esposo and esposa
#> S = 11, p-value = 0.02869
#> alternative hypothesis: true median difference is greater than 0
#> 95 percent confidence interval:
#>    0 Inf
#> sample estimates:
#> median of x-y 
#>             1 
#> 
#> Achieved and Interpolated Confidence Intervals: 
#> 
#>                   Conf.Level L.E.pt U.E.pt
#> Lower Achieved CI     0.9283      0    Inf
#> Interpolated CI       0.9500      0    Inf
#> Upper Achieved CI     0.9755      0    Inf

Los resultados anteriores, también se pueden obtener con la función binom.test de la distribución base de R, como sigue.

binom.test(
  x = s, n = n, p = 0.5, alternative = "greater", 
  conf.level = 0.95
)

#> 
#>  Exact binomial test
#> 
#> data:  s and n
#> number of successes = 11, number of trials = 14, p-value = 0.02869
#> alternative hypothesis: true probability of success is greater than 0.5
#> 95 percent confidence interval:
#>  0.5343434 1.0000000
#> sample estimates:
#> probability of success 
#>              0.7857143

2.2 Prueba de Rangos Asignados de Wilcoxon

2.2.1 Racionalización de la Prueba de Rangos Asignados de Wilcoxon

2.2.2 Implementación en R de la Prueba de Rangos Asignados de Wilcoxon

Ejemplo 2.2 (Prueba de Rangos Asignados de Wilcoxon) Existe considerable evidencia acerca de que los adultos son capaces de utilizar señales visuales en el procesamiento de información auditiva. En una conversación normal, las personas pueden utilizar los movimientos de los labios en el procesamiento de la charla. La congruencia entre los movimientos de los labios y los sonidos del habla son particularmente benéficos en ambientes ruidosos. La investigación ha demostrado que el procesamiento del habla se deteriora cuando las señales auditivas y visuales no son congruentes. En los niños, la habilidad para discriminar y localizar la fuente de estímulos auditivos y visuales complejos se establece alrededor de los seis meses de edad.

Se diseñó un experimento para determinar si los niños de 10 a 16 semanas de edad se dan cuenta de la sincronía entre los movimientos de los labios y los sonidos del habla en una conversación normal. Los niños se colocaron en una habitación a prueba de ruido, que tenía una ventana a través de la cual podían ver a una persona hablando. La persona hablaba en un micrófono y el sonido era dirigido directamente al cuarto (en sincronía) o con una demora de 400 milisegundos (fuera de sincronía). En cada condición se midió el tiempo que el niño miraba la cara de la persona que hablaba. Se argumentó que si el pequeño es capaz de discriminar ambas condiciones, la cantidad de tiempo de ver a la persona sería diferente, aunque a priori no se planteó en cuál de las dos condiciones el tiempo sería mayor (en términos de la poca experiencia en sincronía y lo novedoso que era fuera de sincronía).

Existen considerables diferencias individuales entre los infantes respecto al tiempo que pasaron atendiendo al estímulo. Sin embargo, la diferencia en el tiempo que pasaron viendo en la condición en sincronía y el tiempo que pasaron viendo en la condición fuera de sincronía podría ser un indicador confiable de la capacidad de discriminar. Si el niño pasa más tiempo atendiendo al estímulo en la presentación sincrónica, la diferencia sería positiva; y si el niño pasa más tiempo atendiendo al estímulo en la presentación asincrónica, la diferencia sería negativa. Si el pequeño es capaz de discriminar, las diferencias deberían tender hacia una dirección; más aún, cualquier diferencia en dirección contraria debería ser relativamente pequeña.

Aunque el investigador confía en que las diferencias en el tiempo promedio que se pasaron mirando indican las diferencias en la atención, no está seguro de que las puntuaciones sean suficientemente precisas para que sean representadas en una escala que no sea ordinal. Esto es, sólo puede afirmar que las grandes diferencias reflejan incrementos en la atención; por ejemplo, una diferencia de 30 indica una mayor diferencia en atención que una diferencia de 20. Así, aunque la interpretación de las diferencias en las magnitudes numéricas en el tiempo de mirar no reflejan directamente las diferencias en la atención, el establecer los rangos de las diferencias en el mirar reflejará el orden de las diferencias en el atender al estímulo.

Aquí, la hipótesis nula que se quiere contrastar es que la cantidad de tiempo que pasan los niños viendo a través de la ventana no depende del tipo de presentación. En términos de la prueba de Wilcoxon, la suma de los rangos positivos no difiere de la suma de los rangos negativos. Mientras que, la hipótesis alterna es que la cantidad de tiempo que los niños pasan viendo depende del tipo de presentación; la suma de los rangos positivos difiere de la suma de los rangos negativos. Los datos son tomados de Siegel (1995, pág. 116) y se muestran con el siguiente código.

niños <- tibble(
  sujetos = c(
    "DC", "MK", "VH", "JM", "SB", "MM", "RH", "DJ",
    "JD", "ZC", "CW", "AF"
  ),
  sincronia = c(
    20.3, 17, 6.5, 25, 5.4, 29.2, 2.9, 6.6, 15.8, 8.3,
    34, 8
  ),
  asincronia = c(
    50.4, 87, 25.1, 28.5, 26.9, 36.6, 1, 43.8, 44.2,
    10.4, 29.9, 27.7
  ),
  di = asincronia - sincronia,
  rango = sign(di) * rank(abs(di), ties.method = "average")
)
kableExtra::kbl(
  niños,
  col.names = c(
    "Sujetos", "En sincronía", "Fuera de sincronía",
    "$D_i$", "$R_i$"
  ),
  format.args = list(decimal.mark = ",", big.mark = "."),
  booktabs = TRUE,
  caption = "\\label{tab2:pareja}Porcentaje de falta de atención en presencia de sincronía y sin ella",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  row_spec(which(niños$di == 0), color = "red",
           strikeout = T) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 2.2: Porcentaje de falta de atención en presencia de sincronía y sin ella

Sujetos	En sincronía	Fuera de sincronía	\(D_i\)	\(R_i\)
DC	20,3	50,4	30,1	10
MK	17,0	87,0	70,0	12
VH	6,5	25,1	18,6	6
JM	25,0	28,5	3,5	3
SB	5,4	26,9	21,5	8
MM	29,2	36,6	7,4	5
RH	2,9	1,0	-1,9	-1
DJ	6,6	43,8	37,2	11
JD	15,8	44,2	28,4	9
ZC	8,3	10,4	2,1	2
CW	34,0	29,9	-4,1	-4
AF	8,0	27,7	19,7	7

Una gráfica de caja y bigote, previa a la aplicación de la prueba, podría dar luces respecto a que si los niños son capaces de discriminar entre la sincronía y asincronía. Esto se puede hacer con el siguiente segmento de código.

plot_ly(data = niños, type = "box") %>% 
  add_boxplot(
    y = ~sincronia, name = "En sincronía", jitter = 0.3,
    boxpoints = "all", pointpos = 0,
    marker = list(color = 'rgb(9,56,125)'),
    line = list(color = 'rgb(9,56,125)')
  ) %>% 
  add_boxplot(
    y = ~asincronia, name = "Fuera de sincronía", 
    jitter = 0.3, boxpoints = "all", pointpos = 0,
    marker = list(color = 'rgb(9,56,125)'),
    line = list(color = 'rgb(9,56,125)')
  ) %>%
  layout(
    title = "Tiempo en condición de sincronía y fuera de sincronía",
    showlegend = TRUE
  ) %>%
  layout(
    xaxis = list(title = "Condición", zeroline = F),
    yaxis = list(title = "Tiempo", zeroline = F)
  ) %>%
  config(locale = "es")

El valor del estadístico de la prueba de rangos con signos de Wilcoxon y el p.valor se puede obtener con el siguiente script.

n <- sum(niños$di != 0)
Tmas <- sum(niños$rango[niños$rango > 0])
p.valor <- 2*psignrank(q = Tmas - 1, n = n, lower.tail = FALSE)
paste("El valor del estadístico Tmas =", Tmas)
paste("El p.valor =", round(x = p.valor, digits = 4))

#> [1] "El valor del estadístico Tmas = 73"
#> [1] "El p.valor = 0.0049"

Como se observa en la salida anterior, el \(p_{valor}=73\) es menor que \(\alpha=0,05\), lo que indica que existe suficiente evidencia estadística de que la cantidad de tiempo que los niños pasan viendo a través de la ventana depende del tipo de presentación. Por otro lado, la prueba de rangos con signos de Wilcoxon se puede aplicar, de manera directa, con la función wilcoxon.test del paquete exactRankTests.

library(exactRankTests)
exactRankTests::wilcox.exact(
  x = niños$sincronia, y = niños$asincronia, 
  alternative = "t", paired = T, exact = T, 
  conf.int = T, conf.level = 0.95
)

#> 
#>  Exact Wilcoxon signed rank test
#> 
#> data:  niños$sincronia and niños$asincronia
#> V = 5, p-value = 0.004883
#> alternative hypothesis: true mu is not equal to 0
#> 95 percent confidence interval:
#>  -32.95  -5.45
#> sample estimates:
#> (pseudo)median 
#>        -17.225

La función wilcox.test del paquete stats también realiza esta prueba, como se muestra en el siguiente código.

stats::wilcox.test(
  x = niños$sincronia, y = niños$asincronia,  
  alternative = "t", paired = T, exact = T, 
  conf.int = T, conf.level = 0.95
)

#> 
#>  Wilcoxon signed rank exact test
#> 
#> data:  niños$sincronia and niños$asincronia
#> V = 5, p-value = 0.004883
#> alternative hypothesis: true location shift is not equal to 0
#> 95 percent confidence interval:
#>  -32.95  -5.45
#> sample estimates:
#> (pseudo)median 
#>        -17.225

Ejemplo 2.3 (Prueba de Rangos Asignados de Wilcoxon) Los internos de una prisión federal sirvieron como sujetos en un estudio de toma de decisiones. Primero se midió de manera individual la utilidad de los cigarrillos para los reclusos (valor subjetivo), ya que el cigarrillo es el objeto más negociado en las prisiones. Haciendo uso de la función de utilidad de cada sujeto, el investigador intentó predecir las decisiones que cada hombre haría en un juego donde repetidamente tendría que escoger entre dos opciones y en las cuales pudiera perder o ganar cigarrillos.

La hipótesis evaluada fue que los investigadores podrían predecir las decisiones de los sujetos en términos del valor de la utilidad, en lugar de suponer que la utilidad de los cigarrillos es equivalente al valor nominal de los mismos y, por tanto, predecir la elección “racional” en términos del valor objetivo. Esta hipótesis fue confirmada.

Sin embargo, como era de esperarse, algunas respuestas no se predijeron con éxito, por la hipótesis de maximización de la utilidad esperada. Los investigadores habían conjeturado que tales errores en la predicción se deberían a la probable indiferencia de los sujetos hacia las dos opciones disponibles. Esto es, un recluso podía encontrar igual de atractivo ambas opciones o no parecerle atractiva ninguna y, por lo tanto, ser indiferente a la elección entre ambas opciones. Tales elecciones (de indiferencia) son difíciles de predecir. Pero en tales casos, se razonó que el sujeto vacilaría en plantear una apuesta y tardaría más en decidir. Es decir, la latencia entre el ofrecimiento de las opciones y el aceptar alguna de ellas sería mayor. La segunda hipótesis fue que las latencias o tiempos de respuesta que no se predijeron satisfactoriamente por la maximización de la utilidad esperada, serían mayores que las elecciones predichas con éxito. La hipótesis que se quiere probar aquí es que las latencias de las elecciones predichas incorrectamente serán mayores que las predichas correctamente. Los datos son tomados de Siegel (1995, pág. 119) y se muestran con el siguiente código.

di <- c(
  -2L, 0L, 0L, 1L, 0L, 0L, 4L, 4L, 1L, 1L, 5L, 3L, 5L, 
  3L, -1L, 1L, -1L, 5L, 8L, 2L, 2L, 2L, -3L, -2L, 1L, 
  4L, 8L, 2L, 3L, -1L
)
reclusos <- tibble(
  diferencia = dplyr::if_else(
    di == 0, NA_integer_, di
  )
) %>% 
  mutate(
    rango = sign(diferencia) * 
      rank(
        abs(diferencia), na.last = TRUE, 
             ties.method = "average"
      )
  ) %>% 
  dplyr::transmute(
    reclusos = 1:length(di),
    diferencia = di,
    rango 
  )
options(knitr.kable.NA = "")
kableExtra::kbl(
  reclusos,
  col.names = c(
    "Reclusos", "$D_i$", "$R_i$"
  ),
  format.args = list(decimal.mark = ",", big.mark = "."),
  booktabs = TRUE,
  caption = "\\label{tab2:pareja}Diferencia en el tiempo promedio entre las decisiones correcta e incorrectamente
predichas de los reclusos",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  row_spec(which(reclusos$diferencia == 0), color = "red",
           strikeout = TRUE) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 2.3: Diferencia en el tiempo promedio entre las decisiones correcta e incorrectamente predichas de los reclusos

Reclusos	\(D_i\)	\(R_i\)
1	-2	-11,5
2	0
3	0
4	1	4,5
5	0
6	0
7	4	20,0
8	4	20,0
9	1	4,5
10	1	4,5
11	5	23,0
12	3	16,5
13	5	23,0
14	3	16,5
15	-1	-4,5
16	1	4,5
17	-1	-4,5
18	5	23,0
19	8	25,5
20	2	11,5
21	2	11,5
22	2	11,5
23	-3	-16,5
24	-2	-11,5
25	1	4,5
26	4	20,0
27	8	25,5
28	2	11,5
29	3	16,5
30	-1	-4,5

El siguiente código permite hallar el estadístico y el p.valor de la prueba de Wilcoxon.

Tmas <- reclusos %>% 
  dplyr::select(
    rango
  ) %>% 
  dplyr::filter(
    rango > 0
  ) %>% 
  sum(., na.rm = TRUE)
n <- length(di[di != 0])
p.valor <- psignrank(
  q = Tmas - 1, n = n, lower.tail = FALSE
)
paste("El valor del estadístico Tmas =", Tmas)
paste(
  "El p.valor =", format(
    round(p.valor, digits = 5),
    mode = "double",
    big.mark = ".",
    decimal.mark = ","
  )
)

#> [1] "El valor del estadístico Tmas = 298"
#> [1] "El p.valor = 0,00058"

La salida anterior muestra que el \(p_{valor}=0,000583\), lo cual indica que hay suficiente evidencia estadística para afirmar que las latencias de las elecciones predichas incorrectamente serían mayores que las predichas correctamente.

library(exactRankTests)
wilcox.exact(x = di, y = NULL, alternative = "g", mu = 0,
             paired = FALSE, exact = T)

#> 
#>  Exact Wilcoxon signed rank test
#> 
#> data:  di
#> V = 298, p-value = 0.0005562
#> alternative hypothesis: true mu is greater than 0

En la salida anterior se ha ejecutado la prueba de rangos con signos de Wilcoxon usando la distribución exacta del estadístico \(T^{+}\) de Wilcoxon para muestras relacionadas. A continuación se muestra la ejecución de esta prueba haciendo uso de la convergencia asintótica a la normal de este estadístico.

media <- n*(n + 1)/4
var <- n*(n + 1)*(2*n + 1)/24
z <- (Tmas - media)/sqrt(var)
p.valor.aprox <- pnorm(q = z, lower.tail = F)
paste(
  "El valor del estadístico Z =", round(x = z, digits = 4)
)
paste("El p.valor.aprox =", round(x = p.valor, digits = 4))

#> [1] "El valor del estadístico Z = 3.1113"
#> [1] "El p.valor.aprox = 6e-04"

Como se puede notar el \(p_{valor.exacto} = 0,000583\) , no es muy diferente del \(p_{valor.aprox} = 0,000931\). El resultado anterior puede encontrarse por medio de la función wilcox.exact del paquete exactRankTests, haciendo el argumento exact = F, como sigue.

wilcox.exact(
  x = di, y = NULL, alternative = "g", mu = 0,
  paired = FALSE, exact = F
)

#> 
#>  Asymptotic Wilcoxon signed rank test
#> 
#> data:  di
#> V = 298, p-value = 0.0008779
#> alternative hypothesis: true mu is greater than 0

La salida anterior muestra el \(p_{valor aprox.}\) sin corrección por continuidad. El siguiente código lo muestra usando tal corrección.

z <- (Tmas - 0.5 - media)/sqrt(var)
p.valor.aprox.cor <- pnorm(q = z, lower.tail = F)
paste(
  "El valor del estadístico Z =", round(x = z, digits = 4)
)
paste(
  "El p.valor.aprox.cor =", round(x = p.valor, digits = 4)
)

#> [1] "El valor del estadístico Z = 3.0986"
#> [1] "El p.valor.aprox.cor = 6e-04"

wilcox.test(
  x = di, y = NULL, alternative = "g", mu = 0, 
            paired = F, exact = F, correct = T
)

#> 
#>  Wilcoxon signed rank test with continuity correction
#> 
#> data:  di
#> V = 298, p-value = 0.0009168
#> alternative hypothesis: true location is greater than 0

Note que en los datos de este ejemplo existen rangos empatados, por lo que es necesario ajustar la prueba estadística para considerar el decremento en la variabilidad de \(T^{+}\). La corrección requiere contar los empates y reducir la varianza, respectivamente.

empates <- filter(
  as.data.frame(table(rank(abs(di[di != 0])))), Freq > 1
) %>% 
  dplyr::transmute(
  grupo = 1:nrow(.), 
  rango = Var1,
  Tj = Freq
) %>%
  as.tibble()
kableExtra::kbl(
  empates,
  col.names = c(
    "Grupo", "$R_j$", "$T_j$"
  ),
  format.args = list(decimal.mark = ",", big.mark = "."),
  booktabs = TRUE,
  caption = "\\label{tab2:empates}Grupos de rangos empatados",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2()

Tabla 2.4: Grupos de rangos empatados

Grupo	\(R_j\)	\(T_j\)
1	4.5	8
2	11.5	6
3	16.5	4
4	20	3
5	23	3
6	25.5	2

El cálculo del p.valor corregido por empates se muestra con el siguiente código.

cor.var <- sum(empates$Tj * (empates$Tj^2 - 1)) / 2
var.cor <- var - cor.var
z.cor <- (Tmas - media) / sqrt(var.cor)
p.valor <- pnorm(q = z, lower.tail = F)
paste(
  "El valor del estadístico Z =", round(
    x = z.cor, digits = 6
  )
)
paste("El p.valor =", round(x = p.valor, digits = 6))

#> [1] "El valor del estadístico Z = 3.634119"
#> [1] "El p.valor = 0.000972"

En relación con el cálculo anterior, si se realiza la corrección por continuidad, se obtienen los siguientes resultados.

z.cor.cont <- (Tmas - 0.5 - media) / sqrt(var.cor)
p.valor <- pnorm(q = z, lower.tail = F)
paste("El valor del estadístico Z =", round(
  x = z.cor.cont,
  digits = 6
  )
)
paste("El p.valor =", round(x = p.valor, digits = 6))

#> [1] "El valor del estadístico Z = 3.619286"
#> [1] "El p.valor = 0.000972"

2.3 Prueba para Datos en Forma de Frecuencias de Mcnemar

2.3.1 Racionalización de la Prueba para Datos en Forma de Frecuencias de Mcnemar

2.3.2 Implementación en R de la Prueba para Datos en Forma de Frecuencias de Mcnemar

Ejemplo 2.4 (Prueba para Datos en Forma de Frecuencias de Mcnemar) Durante las campañas presidenciales (y algunas otras campañas para puestos públicos) de 1980 en Estados Unidos se realizaron debates televisivos entre dos o más candidatos. Un investigador en técnicas de comunicación estaba interesado -tanto como los candidatos- en determinar si los debates entre los candidatos presidenciales en las elecciones de 1980 eran efectivos o no en cuanto a cambiar las preferencias de los televidentes hacia los distintos candidatos. Se predijo que si los candidatos (Jimmy Carter y Ronald Reagan) eran igualmente efectivos, habría cambios comparables en las preferencias a cada candidato por parte de los televidentes. Por otro lado, si un candidato era más efectivo o persuasivo durante el debate, entonces habría un cambio diferencial de un candidato a otro. Para evaluar la efectividad del debate, el investigador seleccionó 70 adultos al azar antes del debate y les pidió que indicaran sus preferencias hacia ambos candidatos. Después de la conclusión del debate, les volvió a preguntar acerca de su predilección. Así, en cada caso el conocía las preferencias de las personas antes del debate y después del mismo. Los resultados, obtenidos de Siegel (1995, pág. 103), se presentan en la siguiente tabla. En base a estos datos, pruebe la hipótesis de que entre los televidentes que cambiaron sus preferencias, la probabilidad de que hayan cambiado de Reagan a Carter será la misma de los que cambiaron de Carter a Reagan, con un nivel de confianza del 95%.

carter <- tibble(
  antes = c("Carter", "Reagan"),
  carter = c(28L, 7L),
  regan = c(13L, 27L)
)
kableExtra::kbl(
  carter,
  col.names = c(
    "Preferencia antes del debate televisivo", "Carter",
    "Reagan"
  ),
  align = c("l", "c", "c"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  booktabs = TRUE,
  caption = "\\label{tab2:empates}Preferencias de los sujetos acerca de los candidatos presidenciales antes
y después del debate televisivo",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2() %>%
  column_spec(c(1, 2, 3), width = c("8em", "3em", "3em")) %>% 
  add_header_above(c(" " = 1, "Preferencia después del debate televisivo" = 2))

Tabla 2.5: Preferencias de los sujetos acerca de los candidatos presidenciales antes y después del debate televisivo

	Preferencia después del debate televisivo
Preferencia antes del debate televisivo	Carter	Reagan
Carter	28	13
Reagan	7	27

La prueba exacta de McNemar, usando la distribución binomial, se obtiene con la ejecución del siguiente código.

b <- carter$regan[1]
c <- carter$carter[2]
n <- b + c
p.valor <- 2*pbinom(q = c, size = n, prob = 0.5)
paste("El valor de b =", b)
paste("El valor de c =", c)
paste("El p.valor =", round(x = p.valor, digits = 4))

#> [1] "El valor de b = 13"
#> [1] "El valor de c = 7"
#> [1] "El p.valor = 0.2632"

El cálculo del estadístico \(\chi^{2}\) de Mcnemar y el \(p_{valor}\), usando la convergencia asintótica a la distribución Ji-cuadrada, se muestra con el siguiente código.

X2 <- (abs(b - c) - 1) ^ 2 / (b + c)
p.valor.aprox <- pchisq(q = X2, df = 1, lower.tail = F)
paste("El valor del estadístico X2 =", round(x = X2, digits = 4))
paste("El p.valor.aprox =", round(x = p.valor.aprox, digits = 4))

#> [1] "El valor del estadístico X2 = 1.25"
#> [1] "El p.valor.aprox = 0.2636"

Como se muestra en la salida anterior, dado que el \(p_{valor\,aprox}\) es mayor que el nivel de significancia \(\alpha = 0,05\) no se puede rechazar la hipótesis de que los candidatos fueron igual de efectivos para cambiar la preferencia de los televidentes.

El resultados anterior, se puede obtener de manera directa, con la aplicación de la función mcnemar.test del paquete stats, como sigue.

mcnemar.test(x = as.matrix(carter[ , 2:3]), correct = T)

#> 
#>  McNemar's Chi-squared test with continuity correction
#> 
#> data:  as.matrix(carter[, 2:3])
#> McNemar's chi-squared = 1.25, df = 1, p-value = 0.2636

La función mcnemar.exact del paquete exact2x2, también ejecuta la prueba para datos en forma de frecuencia de Mcnemar. Como se muestra en el siguiente código.

library(exact2x2)
exact2x2::mcnemar.exact(
  x = as.matrix(carter[ , 2:3]), conf.level = 0.95
)

#> 
#>  Exact McNemar test (with central confidence intervals)
#> 
#> data:  as.matrix(carter[, 2:3])
#> b = 13, c = 7, p-value = 0.2632
#> alternative hypothesis: true odds ratio is not equal to 1
#> 95 percent confidence interval:
#>  0.6886514 5.4973372
#> sample estimates:
#> odds ratio 
#>   1.857143

3 Pruebas que Usan Datos de Dos Muestras Independientes

1. Prueba de la Probabilidad Exacta de Fisher* (Sección 3.1)

2. Prueba de la Mediana (Sección 3.2)

3. Prueba \(U\) de Mann-Whitney (Sección ??)

4. Prueba de las Rachas* (Sección ??)

5. Prueba Ji-cuadrado (Sección ??)

6. Prueba de Kolmogorov-Smirnov (Sección ??)

7. Prueba para Rangos Ligados de Moses para Diferencias en la Escala (Sección ??)

3.1 Prueba de la Probabilidad Exacta de Fisher

3.1.1 Racionalización de la Prueba de la Probabilidad Exacta de Fisher

3.1.2 Implementación en R de la Prueba de la Probabilidad Exacta de Fisher

Ejemplo 3.1 (Prueba de la Probabilidad Exacta de Fisher) En un estudio acerca de las situaciones en las cuales las personas amenazan con suicidarse saltando desde un edificio, un puente o una torre, se advirtió que el abucheo o el hostigamiento por parte de una multitud como espectadora ocurría sólo en algunos casos. Varias teorías proponen que un estado psicológico de disminución de la identidad y la autoconciencia, conocido como deindividuación, puede contribuir al fenómeno de hostigamiento. Se conocen algunos factores que inducen reacciones en las multitudes, incluidos la temperatura, el ruido y la fatiga. En un esfuerzo por evaluar varias hipótesis concernientes al hostigamiento por parte de las multitudes, Mann revisó 21 artículos publicados acerca de suicidio y examinó la relación entre el hostigamiento por parte de la multitud y el mes del año; esto último se refería más bien al índice de temperatura. La hipótesis es que habría un incremento en el hostigamiento por parte de los espectadores cuando hiciera calor. Prueba la hipótesis anterior con un nivel de significancia del 10%. Los datos fueron tomados de Siegel (1995, pág. 137). El siguiente código de R, muestra la tabla de datos.

datos4.1 <- tibble(
  mes = c(
    "Junio-septiembre", "Octubre-mayo", "Total"
  ),
  hostigamiento = c(8L, 2L, 10L),
  no_hostigamiento = c(4L, 7L, 11L),
  combinacion = hostigamiento + no_hostigamiento
)
kableExtra::kbl(
  datos4.1,
  col.names = c(
    "Mes", "Hostigamiento", "No hostigamientp",
    "Combinación"
  ),
  align = c("l", "c", "c", "c"),
  format.args = list(decimal.mark = ",", big.mark = "."),
  booktabs = TRUE,
  caption = "\\label{tab2:hostigamiento}Incidencia de hostigamiento en episodios de intento de suicidio",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = NULL,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2() %>%
  column_spec(c(1, 2, 3, 4), width = c("8em", "5em", "5em",
                                    "5em")) %>% 
  add_header_above(c(" " = 1, "Multitud" = 2, " " = 1)) %>% 
  column_spec(c(1, 3), border_right = c(TRUE, TRUE))

Tabla 3.1: Incidencia de hostigamiento en episodios de intento de suicidio

	Multitud
Mes	Hostigamiento	No hostigamientp	Combinación
Junio-septiembre	8	4	12
Octubre-mayo	2	7	9
Total	10	11	21

El ||{valor}$ de la prueba de la probabilidad exacta de Fisher, se obtiene a partir del siguiente segmento de código. En esta salida se puede observar que el \(p_{valor}\) (0.0563) es menor que el nivel de significancia especificado (10%), lo que indica que el hostigamiento por parte de las multitudes que atestiguan una amenaza de suicidio es afectada por la temperatura (medida según el mes del año).

p.valor <- phyper(
  q = as.integer(datos4.1[1, 2] - 1), 
  m = as.integer(datos4.1[3, 2]),
  n = as.integer(datos4.1[3, 3]), 
  k = as.integer(datos4.1[1, 4]),
  lower.tail = FALSE
)
paste(
  "El p.valor =", 
  format(
    round(p.valor, digits = 6),
    mode = "double",
    big.mark = ".",
    decimal.mark = ","
  )
)

#> [1] "El p.valor = 0,056323"

La función fisher.test del paquete stats, permite ejecutar, de manera directa, la prueba de la probabilidad exacta de Fisher, como se muest/14en el siguiente código. Fíjese que el \(p_{valor}\) obtenido en esta salida es idéntico al anterior.

fisher.test(x = datos4.1[1:2, 2:3], 
            alternative = "greater",
            conf.int = TRUE , conf.level = 0.9)

#> 
#>  Fisher's Exact Test for Count Data
#> 
#> data:  datos4.1[1:2, 2:3]
#> p-value = 0.05632
#> alternative hypothesis: true odds ratio is greater than 1
#> 90 percent confidence interval:
#>  1.293097      Inf
#> sample estimates:
#> odds ratio 
#>   6.302622

Otra función que permite ejecutar esta prueba es fisher.exact, la cual se encuentra en el paquete exact2x2. El siguiente código muestra la aplicación de la prueba por medio de la mencionada función.

library(exact2x2)
fisher.exact(x = datos4.1[1:2, 2:3], 
             alternative = "greater",
             conf.int = TRUE, conf.level = 0.9)

#> 
#>  One-sided Fisher's Exact Test
#> 
#> data:  datos4.1[1:2, 2:3]
#> p-value = 0.05632
#> alternative hypothesis: true odds ratio is greater than 1
#> 90 percent confidence interval:
#>  1.293097      Inf
#> sample estimates:
#> odds ratio 
#>   6.302622

El \(p_{valor}\) aproximado se puede obtener usando la aproximación asintótica a la distribución Ji-cuadrada, con el siguiente script.

X2 <- (as.integer(datos4.1[3, 4]) * (abs(as.integer(datos4.1[1, 2]) * as.integer(datos4.1[2, 3]) - as.integer(datos4.1[1, 3]) * as.integer(datos4.1[2, 2])) - as.integer(datos4.1[3, 4]) / 2)^2) / (as.integer(datos4.1[1, 4]) * as.integer(datos4.1[2, 4]) * as.integer(datos4.1[3, 2]) * as.integer(datos4.1[3, 3]))
p.valor <- pchisq(q = X2, df = 1, lower.tail = FALSE) / 2
paste(
  "El p.valor =", 
  format(
    round(p.valor, digits = 6),
    mode = "double",
    big.mark = ".",
    decimal.mark = ","
  )
)

#> [1] "El p.valor = 0,057439"

3.2 Prueba de la Mediana

3.2.1 Racionalización de la Prueba de la Mediana

3.2.2 Implementación en R de la Prueba de la Mediana

Ejemplo 3.2 (Prueba de la Mediana) Los datos siguientes representan la duración en horas de dos tipos diferentes de baterías: \[ \begin{array}{ccccccc} \hline \mathbf{Tipo \, I}: 40.1 & 30.5 & 40.2 & 45 & 55.5 & 30.6\\ \mathbf{Tipo \, II}: 50.5 & 50.7 & 45.5 & 55.7 & 60 & 40.8\\ \hline \end{array} \] Contraste la hipótesis, con un nivel de riesgo del 5%, de que la mediana de la duración de la batería del tipo I es menor que la del tipo II. Los datos fueron tomados de Gómez (2005, pág. 437).

Para ejecutar la prueba se deben ordenar las muestras combinadas de manera creciente para después identificar el número observaciones que son menores que la mediana global \(\left ( M = 45,25 \right )\) y las que son mayores o iguales que este valor para cada tipo de batería, como se muestra a continuación, con el siguiente código.

I <- c(40.1, 30.5, 40.2, 45, 55.5, 30.6)
II <- c(50.5, 50.7, 45.5, 55.7, 60, 40.8)
Dur <- c(I, II)
TB <- factor(rep(x = c("I", "II"), 
                 times = c(length(I), length(II))))
M <- median(Dur)
baterias <- tibble(
  duracion = Dur,
  tb = TB,
  comparacion = factor(
    if_else(Dur < M, "$< M$", "$\\ge M$")
  )
) %>% 
  arrange(duracion)
kableExtra::kbl(
  baterias,
  col.names = c(
    "Duración", "Tipo de batería", "Comparación"
  ),
  format.args = list(decimal.mark = ",", big.mark = "."),
  booktabs = TRUE,
  caption = "\\label{tab2:baterias}Tiempo de duración de las baterías (horas)",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = FALSE,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2() %>%
  scroll_box(width = "100%", height = "400px")

Tabla 3.2: Tiempo de duración de las baterías (horas)

Duración	Tipo de batería	Comparación
30,5	I	\(< M\)
30,6	I	\(< M\)
40,1	I	\(< M\)
40,2	I	\(< M\)
40,8	II	\(< M\)
45,0	I	\(< M\)
45,5	II	\(\ge M\)
50,5	II	\(\ge M\)
50,7	II	\(\ge M\)
55,5	I	\(\ge M\)
55,7	II	\(\ge M\)
60,0	II	\(\ge M\)

Ahora, el siguiente código muestra la tabla de contingencia con el número observaciones que son menores y mayores o iguales que la mediana global para cada tipo de batería.

tf_baterias <- baterias %>% 
  group_by(tb, comparacion) %>% 
  tally() %>% 
  spread(tb, n) %>% 
  dplyr::arrange(desc(comparacion)) %>% 
  dplyr::mutate(
    Total = I + II
  ) %>% 
  dplyr::select_if(is.numeric) %>% 
  add_row(., summarise_all(., sum)) %>% 
  dplyr::mutate(
    comparacion = c("$ < M $", "$ \\ge M $", "Total")
  ) %>% 
  dplyr::select(comparacion, I, II, Total)
kableExtra::kbl(
  tf_baterias,
  col.names = c(
    "Comparación", "Tipo I", "Tipo II", "Total"
  ),
  format.args = list(decimal.mark = ",", big.mark = "."),
  booktabs = TRUE,
  caption = "\\label{tab2:baterias}Duración de baterias menores y mayores o iguales que la mediana global",
  escape = FALSE
) %>%
  kable_styling(
    bootstrap_options = "striped",
    full_width = FALSE,
    fixed_thead = TRUE
  ) %>%
  kable_classic_2("hover") %>% 
  add_header_above(c(" " = 1, "Baterías" = 2, " " = 1)) %>% 
  column_spec(c(1, 3), border_right = c(TRUE, TRUE))  

Tabla 3.3: Duración de baterias menores y mayores o iguales que la mediana global

	Baterías
Comparación	Tipo I	Tipo II	Total
$ < M $	5	1	6
$ M $	1	5	6
Total	6	6	12

En este orden de ideas, el siguiente procedimiento muestra el estadístico y el p.valor de la prueba de la mediana de Mood-Westenberg.

p.valor <- phyper(
  q = as.integer(tf_baterias[1, 2]) - 1, 
  m = as.integer(tf_baterias[3, 2]), 
  n = as.integer(tf_baterias[3, 3]),
  k = as.integer(tf_baterias[1, 4]), 
  lower.tail = FALSE
)
paste("El valor del Estadístico U =", tf_baterias[1, 2])
paste("El p.valor =", round(x = p.valor, digits = 5))

#> [1] "El valor del Estadístico U = 5"
#> [1] "El p.valor = 0.04004"

De la salida anterior se observa que el \(p.valor =0,040043\), lo cual indica que las batería del tipo I duran menos que las del tipo II.

library(RVAideMemoire)
mood.medtest(x = Dur, g = TB, exact = T)

#> 
#>  Mood's median test
#> 
#> data:  Dur by TB
#> p-value = 0.08009

La prueba exacta de Fisher se puede ejecutar con la función fisher.test, como se muestra en la siguiente salida, generada con el siguiente código.

fisher.test(x =  select_if(tf_baterias[1:2, 2:3], is.numeric), 
            alternative = "greater", conf.int = 0.95)

#> 
#>  Fisher's Exact Test for Count Data
#> 
#> data:  select_if(tf_baterias[1:2, 2:3], is.numeric)
#> p-value = 0.04004
#> alternative hypothesis: true odds ratio is greater than 1
#> 95 percent confidence interval:
#>  1.118505      Inf
#> sample estimates:
#> odds ratio 
#>   16.60315

Las funciones mood.medtest del paquete RVAideMemoire muestra el \(p_{valor}\) de la prueba de la mediana de Mood para dos colas. De tal manera que si éste \(p_{valor}\) se divide por dos, se obtiene el p.valor para la prueba a una cola.

RVAideMemoire::mood.medtest(x = Dur, g = TB, exact = T)

#> 
#>  Mood's median test
#> 
#> data:  Dur by TB
#> p-value = 0.08009

El resultado anterior, también se puede obtener con la función mood.median.test del paquete MultNonParam.

MultNonParam::mood.median.test(x = I, y = II, exact = T)

#> $p.value
#> [1] 0.08008658

3.3 Información de Sesión

as_tibble(devtools::session_info()$packages) %>%
  dplyr::select(
    package, loadedversion, source
  ) %>%
  dplyr::filter(
    package %in% packages
  ) %>%
  DT::datatable(
    rownames = FALSE,
    colnames = c("Paquete", "Versión", "Fuente"),
    caption = "Información de sesión",
    filter = "top",
    selection = "multiple",
    class = "cell-border stripe"
  )

Bibliografía

[1] J. Allaire, Y. Xie, J. McPherson, et al. rmarkdown: Dynamic Documents for R. R package version 2.14. 2022. <URL: https://CRAN.R-project.org/package=rmarkdown>.

[2] A. T. Arnholt and B. Evans. BSDA: Basic Statistics and Data Analysis. R package version 1.2.1. 2021. <URL: https://CRAN.R-project.org/package=BSDA>.

[3] S. M. Bache and H. Wickham. magrittr: A Forward-Pipe Operator for R. R package version 2.0.1. 2020. <URL: https://CRAN.R-project.org/package=magrittr>.

[4] C. Boettiger. knitcitations: Citations for Knitr Markdown Files. R package version 1.0.12. 2021. <URL: https://github.com/cboettig/knitcitations>.

[5] G. R. W. a. Bolker, T. Lumley, R. C. Johnson, et al. gmodels: Various R Programming Tools for Model Fitting. R package version 2.18.1. 2018. <URL: https://CRAN.R-project.org/package=gmodels>.

[6] F. Caeiro and A. Mateus. randtests: Testing randomness in R. R package version 1.0. 2014. <URL: https://CRAN.R-project.org/package=randtests>.

[7] L. Carvalho. “An Improved Evaluation of Kolmogorov’s Distribution”. In: Journal of Statistical Software 65.3 (2015), pp. 1-7. <URL: http://www.jstatsoft.org/v65/c03/>.

[8] L. Carvalho. kolmim: An Improved Evaluation of Kolmogorov’s Distribution. R package version 1.0. 2015. <URL: https://CRAN.R-project.org/package=kolmim>.

[9] S. Chasalow. combinat: combinatorics utilities. R package version 0.0-8. 2012. <URL: https://CRAN.R-project.org/package=combinat>.

[10] J. Cheng, C. Sievert, B. Schloerke, et al. htmltools: Tools for HTML. R package version 0.5.2. 2021. <URL: https://github.com/rstudio/htmltools>.

[11] M. P. Fay. “Confidence intervals that match Fisher’s exact or Blaker’s exact tests”. In: Biostatistics 11.2 (2010), pp. 373-374. <URL: https://www.niaid.nih.gov/about/brb-staff-fay>.

[12] M. P. Fay. exactci: Exact P-Values and Matching Confidence Intervals for Simple Discrete Parametric Cases. R package version 1.4-2. 2021. <URL: https://CRAN.R-project.org/package=exactci>.

[13] M. P. Fay. ssanv: Sample Size Adjusted for Nonadherence or Variability of Input Parameters. R package version 1.1. 2015. <URL: https://CRAN.R-project.org/package=ssanv>.

[14] M. P. Fay. “Two-sided Exact Tests and Matching Confidence Intervals for Discrete Data”. In: R Journal 2.1 (2010), pp. 53-58. <URL: https://journal.r-project.org/>.

[15] M. P. Fay. “Two-sided Exact Tests and Matching Confidence Intervals for Discrete Data”. In: R Journal 2.1 (2010), pp. 53-58. <URL: https://journal.R-project.org/>.

[16] M. P. Fay, M. E. Halloran, and D. A. Follmann. “Accounting for variability in sample size estimation with applications to nonadherence and estimation of variance and effect size”. In: Biometrics 63 (2007), pp. 465-474.

[17] M. P. Fay, S. A. Hunsberger, M. Nason, et al. exact2x2: Exact Tests and Confidence Intervals for 2x2 Tables. R package version 1.6.5. 2020. <URL: https://CRAN.R-project.org/package=exact2x2>.

[18] M. P. Fay and K. Lumbard. “Confidence Intervals for Difference in Proportions for Matched Pairs Compatible with Exact McNemars or Sign Tests”. In: unpublished manuscript (2020).

[19] M. P. Fay, M. A. Proschan, and E. Brittain. “Combining one sample confidence procedures for inference in the two sample case”. In: Biometrics (2014). <URL: https://www.niaid.nih.gov/about/brb-staff-fay>.

[20] E. E. Gabriel, M. Nason, M. P. Fay, et al. “A boundary-optimized rejection region test for the two-sample binomial problem”. In: Statistics in Medicine 37.7 (2018), pp. 1047-1058.

[21] A. Genz and F. Bretz. Computation of Multivariate Normal and t Probabilities. Lecture Notes in Statistics. Heidelberg: Springer-Verlag, 2009. ISBN: 978-3-642-01688-2.

[22] A. Genz, F. Bretz, T. Miwa, et al. mvtnorm: Multivariate Normal and t Distributions. R package version 1.1-1. 2020. <URL: http://mvtnorm.R-forge.R-project.org>.

[23] R. K. S. Hankin. “Additive integer partitions in R”. In: Journal of Statistical Software, Code Snippets 16 (1 may. 2006).

[24] R. K. S. Hankin. partitions: Additive Partitions of Integers. R package version 1.10-2. 2021. <URL: https://github.com/RobinHankin/partitions>.

[25] F. E. Harrell Jr. Hmisc: Harrell Miscellaneous. R package version 4.5-0. 2021. <URL: https://CRAN.R-project.org/package=Hmisc>.

[26] L. Henry and H. Wickham. purrr: Functional Programming Tools. R package version 0.3.4. 2020. <URL: https://CRAN.R-project.org/package=purrr>.

[27] M. Hervé. RVAideMemoire: Testing and Plotting Procedures for Biostatistics. R package version 0.9-80. 2021. <URL: https://CRAN.R-project.org/package=RVAideMemoire>.

[28] T. Hothorn and K. Hornik. exactRankTests: Exact Distributions for Rank and Permutation Tests. R package version 0.8-35. 2022. <URL: https://CRAN.R-project.org/package=exactRankTests>.

[29] S. Jankowski. MultNonParam: Multivariate Nonparametric Methods. R package version 1.3.6. 2021. <URL: https://CRAN.R-project.org/package=MultNonParam>.

[30] K. Müller and H. Wickham. tibble: Simple Data Frames. R package version 3.1.5. 2021. <URL: https://CRAN.R-project.org/package=tibble>.

[31] K. Nordhausen, H. Oja, and D. E. Tyler. ICS: Tools for Exploring Multivariate Data via ICS/ICA. R package version 1.3-1. 2018. <URL: https://CRAN.R-project.org/package=ICS>.

[32] K. Nordhausen, H. Oja, and D. E. Tyler. “Tools for Exploring Multivariate Data: The Package ICS”. In: Journal of Statistical Software 28.6 (2008), pp. 1-31. <URL: http://www.jstatsoft.org/v28/i06/>.

[33] K. Nordhausen, S. Sirkia, H. Oja, et al. ICSNP: Tools for Multivariate Nonparametrics. R package version 1.1-1. 2018. <URL: https://CRAN.R-project.org/package=ICSNP>.

[34] P. O. Perry. utf8: Unicode Text Processing. R package version 1.2.1. 2021. <URL: https://CRAN.R-project.org/package=utf8>.

[35] D. Qiu. snpar: Supplementary Non-parametric Statistics Methods. R package version 1.0. 2014. <URL: https://CRAN.R-project.org/package=snpar>.

[36] R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria, 2020. <URL: https://www.R-project.org/>.

[37] B. Ripley. MASS: Support Functions and Datasets for Venables and Ripley’s MASS. R package version 7.3-54. 2021. <URL: http://www.stats.ox.ac.uk/pub/MASS4/>.

[38] B. Rudis. hrbrthemes: Additional Themes, Theme Components and Utilities for ggplot2. R package version 0.8.0. 2020. <URL: http://github.com/hrbrmstr/hrbrthemes>.

[39] D. Sarkar. Lattice: Multivariate Data Visualization with R. ISBN 978-0-387-75968-5. New York: Springer, 2008. <URL: http://lmdvr.r-forge.r-project.org>.

[40] D. Sarkar. lattice: Trellis Graphics for R. R package version 0.20-44. 2021. <URL: http://lattice.r-forge.r-project.org/>.

[41] G. Schneider, E. Chicken, and R. Becvarik. NSM3: Functions and Datasets to Accompany Hollander, Wolfe, and Chicken - Nonparametric Statistical Methods, Third Edition. R package version 1.16. 2021. <URL: https://CRAN.R-project.org/package=NSM3>.

[42] C. Sievert. Interactive Web-Based Data Visualization with R, plotly, and shiny. Chapman and Hall/CRC, 2020. ISBN: 9781138331457. <URL: https://plotly-r.com>.

[43] C. Sievert, C. Parmer, T. Hocking, et al. plotly: Create Interactive Web Graphics via plotly.js. R package version 4.10.0. 2021. <URL: https://CRAN.R-project.org/package=plotly>.

[44] Terry M. Therneau and Patricia M. Grambsch. Modeling Survival Data: Extending the Cox Model. New York: Springer, 2000. ISBN: 0-387-98784-3.

[45] T. M. Therneau. survival: Survival Analysis. R package version 3.2-11. 2021. <URL: https://github.com/therneau/survival>.

[46] R. Vaidyanathan, Y. Xie, J. Allaire, et al. htmlwidgets: HTML Widgets for R. R package version 1.5.4. 2021. <URL: https://github.com/ramnathv/htmlwidgets>.

[47] W. N. Venables and B. D. Ripley. Modern Applied Statistics with S. Fourth. ISBN 0-387-95457-0. New York: Springer, 2002. <URL: https://www.stats.ox.ac.uk/pub/MASS4/>.

[48] H. Wickham. forcats: Tools for Working with Categorical Variables (Factors). R package version 0.5.1. 2021. <URL: https://CRAN.R-project.org/package=forcats>.

[49] H. Wickham. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York, 2016. ISBN: 978-3-319-24277-4. <URL: https://ggplot2.tidyverse.org>.

[50] H. Wickham. reshape2: Flexibly Reshape Data: A Reboot of the Reshape Package. R package version 1.4.4. 2020. <URL: https://github.com/hadley/reshape>.

[51] H. Wickham. “Reshaping Data with the reshape Package”. In: Journal of Statistical Software 21.12 (2007), pp. 1-20. <URL: http://www.jstatsoft.org/v21/i12/>.

[52] H. Wickham. stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.4.0. 2019. <URL: https://CRAN.R-project.org/package=stringr>.

[53] H. Wickham. “testthat: Get Started with Testing”. In: The R Journal 3 (2011), pp. 5-10. <URL: https://journal.r-project.org/archive/2011-1/RJournal_2011-1_Wickham.pdf>.

[54] H. Wickham. testthat: Unit Testing for R. R package version 3.0.4. 2021. <URL: https://CRAN.R-project.org/package=testthat>.

[55] H. Wickham. tidyr: Tidy Messy Data. R package version 1.1.4. 2021. <URL: https://CRAN.R-project.org/package=tidyr>.

[56] H. Wickham. tidyverse: Easily Install and Load the Tidyverse. R package version 1.3.1. 2021. <URL: https://CRAN.R-project.org/package=tidyverse>.

[57] H. Wickham, M. Averick, J. Bryan, et al. “Welcome to the tidyverse”. In: Journal of Open Source Software 4.43 (2019), p. 1686. DOI: 10.21105/joss.01686.

[58] H. Wickham, J. Bryan, and M. Barrett. usethis: Automate Package and Project Setup. R package version 2.1.3. 2021. <URL: https://CRAN.R-project.org/package=usethis>.

[59] H. Wickham, W. Chang, L. Henry, et al. ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. R package version 3.3.5. 2021. <URL: https://CRAN.R-project.org/package=ggplot2>.

[60] H. Wickham, R. François, L. Henry, et al. dplyr: A Grammar of Data Manipulation. R package version 1.0.6. 2021. <URL: https://CRAN.R-project.org/package=dplyr>.

[61] H. Wickham and J. Hester. readr: Read Rectangular Text Data. R package version 1.4.0. 2020. <URL: https://CRAN.R-project.org/package=readr>.

[62] H. Wickham, J. Hester, and W. Chang. devtools: Tools to Make Developing R Packages Easier. R package version 2.4.2. 2021. <URL: https://CRAN.R-project.org/package=devtools>.

[63] Y. Xie. bookdown: Authoring Books and Technical Documents with R Markdown. ISBN 978-1138700109. Boca Raton, Florida: Chapman and Hall/CRC, 2016. <URL: https://bookdown.org/yihui/bookdown>.

[64] Y. Xie. bookdown: Authoring Books and Technical Documents with R Markdown. R package version 0.26. 2022. <URL: https://CRAN.R-project.org/package=bookdown>.

[65] Y. Xie, J. Allaire, and G. Grolemund. R Markdown: The Definitive Guide. ISBN 9781138359338. Boca Raton, Florida: Chapman and Hall/CRC, 2018. <URL: https://bookdown.org/yihui/rmarkdown>.

[66] Y. Xie, J. Cheng, and X. Tan. DT: A Wrapper of the JavaScript Library DataTables. R package version 0.19. 2021. <URL: https://github.com/rstudio/DT>.

[67] Y. Xie, C. Dervieux, and E. Riederer. R Markdown Cookbook. ISBN 9780367563837. Boca Raton, Florida: Chapman and Hall/CRC, 2020. <URL: https://bookdown.org/yihui/rmarkdown-cookbook>.

[68] A. Zeileis and Y. Croissant. “Extended Model Formulas in R: Multiple Parts and Multiple Responses”. In: Journal of Statistical Software 34.1 (2010), pp. 1-13. DOI: 10.18637/jss.v034.i01.

[69] A. Zeileis and Y. Croissant. Formula: Extended Model Formulas. R package version 1.2-4. 2020. <URL: https://CRAN.R-project.org/package=Formula>.

[70] P. Zhao. bookdownplus: Generate Assorted Books and Documents with R bookdown Package. R package version 1.5.8. 2020. <URL: https://github.com/pzhaonet/bookdownplus>.

[71] H. Zhu. kableExtra: Construct Complex Table with kable and Pipe Syntax. R package version 1.3.4. 2021. <URL: https://CRAN.R-project.org/package=kableExtra>.